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Abstract 

We  develop  a  new  continuum  mechanics  modeling  framework  for  liquid-vapor  flows,  with 
particular  focus  on  the  van  der  Waals  fluid.  By  invoking  microforce  theory,  the  Coleman-Noll 
procedure  is  generalized  to  derive  consistent  constitutive  relations  in  the  presence  of  non¬ 
local  effects.  A  new  thermodynamically  consistent  algorithm  for  the  van  der  Waals  model  is 
designed,  using  functional  entropy  variables  and  a  new  temporal  scheme  employing  a  family 
of  new  quadrature  rules.  We  show  that  the  resulting  fully  discrete  scheme  is  unconditionally 
stable  in  entropy  and  second-order  time-accurate.  Isogeometric  analysis  is  utilized  for  spatial 
discretization.  The  analytical  properties  of  the  formulation  are  corroborated  by  benchmark 
problems.  Three  sets  of  application  problems  are  simulated  to  demonstrate  the  capability 
of  the  model  and  the  algorithm.  Our  methodology  provides  a  particularly  useful  predictive 
tool  for  boiling  flows. 

Keywords:  Phase-field  model,  Diffuse  interface,  Microforce,  Coleman-Noll  approach,  Van 
der  Waals  fluid,  Non-convex  flux,  Entropy  variables,  Time  integration,  Isogeometric  analysis, 
Phase  transition,  Evaporation,  Condensation,  Boiling 
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1  Introduction 


1.1  Phase  transition  and  classical  modeling  techniques 

Liquid-vapor  two-phase  flows  are  ubiquitous  in  the  natural  world  as  well  as  in  industry. 
Liquid-vapor  phase  transitions  involve  a  sharp  change  of  the  fluid  density.  Typical  environ¬ 
mental  changes  include  pressure  variations  and  thermal  variations.  For  instance,  the  local 
pressure  near  a  rotating  propeller  may  drop  below  the  boiling  pressure,  and  vapor  bubbles 
may  generate  near  blades  [40].  This  is  the  phenomenon  of  cavitation  and  it  is  still  a  limiting 
factor  for  ship  propeller  design.  Phase  transitions  induced  by  temperature  variations  can  be 
observed  in  daily  life  as  boiling,  evaporation,  and  condensation.  In  industry,  liquid-vapor 
phase  transitions  take  place  in  steam  generators,  heat  exchangers,  and  various  pipelines. 
The  accompanying  thermal  effects  make  multiphase  flow  a  widely-used  mechanism  for  en¬ 
ergy  transfer. 

To  date,  many  modeling  techniques  have  been  designed  to  simulate  multiphase  flows. 
Most  of  them  fall  into  the  categories  of  interface-tracking  methods  or  interface-capturing 
methods.  Interface-tracking  methods  resolve  the  interface  by  aligning  the  computational 
mesh  along  the  interface  and  by  updating  the  mesh  accordingly  with  the  fluid  flow.  This 
approach  gives  a  sharp  and  accurate  representation  of  the  interface.  However,  it  requires  con¬ 
stant  re-meshing  of  the  computational  domain,  and  it  is  typically  intractable  for  topological 
transitions.  Three-dimensional  problems  with  severe  topological  transitions  are  notoriously 
difficult  to  solve  with  interface-tracking  methods.  Interface-capturing  methods  use  additional 
unknowns  to  implicitly  represent  the  interface.  The  interface  is  typically  immersed  in  the 
computational  domain.  Consequently,  the  interface  representation  is  less  accurate  than  that 
of  the  interface-tracking  methods.  However,  the  interface-capturing  methods  enjoy  several 
advantages:  they  are  relatively  easy  to  implement,  the  mesh  updating  burden  is  reduced,  and 
topological  transitions  are  easily  handled.  Existing  instantiations  of  the  interface-capturing 
methods  include  the  volume-of-fluid  (VOF)  method  [31]  and  the  level-set  method  [55].  Both 
methods  have  been  utilized  in  commercial  codes  and  remain  popular  subjects  in  the  litera¬ 
ture.  However,  they  are  not  without  shortcomings.  The  VOF  method  uses  a  post-processing 
procedure  to  construct  the  interface,  which  inevitably  introduces  errors.  In  level-set  methods, 
the  level-set  function  needs  to  be  reinitialized  every  few  steps.  The  reinitialization  procedure 
is  rather  ad  hoc  and  does  not  maintain  the  conservation  structure  of  the  governing  equations. 
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1.2  Phase-field  models 

To  address  the  aforementioned  modeling  difficulties,  phase-field  models  were  proposed  as  an 
alternative  interface-capturing  method.  The  phase-held  method  uses  an  order  parameter  to 
distinguish  different  phases  [1].  Phase-held  models  postulate  that  the  interface  has  finite 
width  and  material  properties  transit  across  the  interfacial  region  smoothly  but  sharply. 
Based  on  these  postulates,  van  der  Waals  developed  his  Nobel  Prize-winning  theory  to  cal¬ 
culate  the  capillarity  for  liquid- vapor  interfaces  [58].  Later,  Korteweg  developed  the  so-called 
Korteweg  stress  formulation  and  coupled  the  van  der  Waals  theory  with  hydrodynamics  [39] . 
The  fluid  equations  based  on  the  van  der  Waals  theory  are  called  the  Navier- Stokes- Korteweg 
equations.  This  theory  is  characterized  by  a  non-convex  free  energy,  supplemented  with  a 
non-local  density  gradient-squared  term.  This  non-local  term  regularizes  the  singularity  in¬ 
troduced  by  the  non-convex  free  energy  function.  The  non-local  term  represents  the  surface 
energy.  In  contemporary  continuum  mechanics,  this  model  falls  into  the  category  of  the 
grade- N  fluid  model  [70].  In  1985,  Dunn  and  Serrin  studied  the  thermodynamic  consistency 
of  the  Navier-Stokes-Korteweg  equations  and  found  that  for  the  model  to  be  consistent  with 
the  second  law,  a  new  term  had  to  be  added  to  the  energy  equation  [17].  They  called  this 
term  the  “interstitial  working  flux” . 

Parallel  to  the  development  of  the  Navier-Stokes-Korteweg  equations,  another  branch  of 
phase-held  models  has  been  developed  focusing  on  multicomponent  systems  [1].  The  original 
idea  emanates  from  the  work  of  Calm  and  Hilliard  [9],  in  which  a  fourth-order  nonlinear 
diffusion  equation  is  proposed  to  mimic  the  behavior  of  a  two-component  mixture.  Recently, 
the  Cahn-Hilliard  type  models  have  been  generalized  to  more  complicated  multicomponent 
systems,  such  as  spinodal  decomposition  [44],  tumor  growth  [72],  fingering  effect  in  porous 
medium  [23],  topology  optimization  [13],  etc.  Significant  progress  was  made  by  Gurtin 
and  his  collaborators  in  providing  a  rational  mechanics  framework  for  the  Cahn-Hilliard 
type  models  [27].  In  Gurtin’s  theory,  microforces  were  introduced  to  account  for  the  phase 
dynamics.  Later,  this  theory  was  applied  to  construct  a  plasticity  theory  of  single  crystals 
[28],  fracture  models  [74],  alloy  models  [47],  and  ferroelectric  models  [69],  to  list  a  few.  In 
Section  2,  this  theory  is  adopted  as  a  means  to  derive  constitutive  relations  for  the  van  der 
Waals  fluid  material.  Interestingly,  the  “interstitial  working  flux”  appears  naturally  in  this 
derivation  as  the  power  expenditure  of  the  microstress.  This,  in  part,  justifies  the  work  of 
Dunn  and  Serrin  within  the  rational  thermo  mechanics  framework  of  [11]. 

Traditional  interface-tracking  and  interface-capturing  methods  are  designed  based  on  geo¬ 
metrical  information  of  existing  interfaces.  When  dealing  with  phase  transition  phenomena, 
those  methods  become  intractable.  One  may  need  to  introduce  artificial  procedures  and 
empirical  assumptions  [68]  to  mimic  such  phenomena.  In  contrast,  the  solid  mathematical 
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and  thermodynamic  foundations  of  phase-field  models  allow  them  to  describe  these  compli¬ 
cated  phenomena  without  resorting  to  modeling  tricks.  In  this  work  this  advantage  will  be 
demonstrated  by  a  suite  of  boiling  simulations.  Boiling  is  regarded  to  be  highly  difficult  for 
numerical  simulations  [3,  16,  38,  53,  41].  Traditional  models  require  empirical  knowledge, 
such  as  bubble  release  rates  and  bubble  departure  radii,  to  describe  the  boiling  process.  In 
this  work,  two-  and  three-dimensional  boiling  simulations  are  carried  out  using  the  van  der 
Waals  fluid  model.  Owing  to  the  thermodynamically  consistent  nature,  the  dependency  on 
empirical  knowledge  is  significantly  reduced,  and  there  are  no  ad  hoc  procedures  involved. 
Our  approach  provides  a  unified  predictive  capability  for  both  nucleate  and  him  boiling 
Despite  their  success,  phase-held  models  face  several  challenges.  The  entropy  function 
for  phase-held  models  are  always  non-convex,  which  creates  difficulty  for  both  mathematical 
analysis  and  numerical  simulation.  Phase-held  models  usually  have  a  high-order  differential 
term,  which  necessitates  novel  numerical  discretization  techniques,  and  the  interface  width 
for  real  materials  is  typically  a  few  nanometers.  Therefore,  adaptive  refinement  near  the 
interfacial  region  is  required  for  real-world  simulations. 

1.3  Numerical  analysis 

In  the  numerical  analysis  of  nonlinear  problems,  a  central  issue  is  stability.  One  relevant 
example  for  the  current  work  is  the  study  of  entropy-stable  schemes  for  gas  dynamics.  It  was 
revealed  that  the  weak  form  of  the  compressible  Navier-Stokes  equations  will  automatically 
satisfy  the  Clausius-Duhem  inequality  by  utilizing  entropy  variables  [34],  In  the  late  1980s, 
the  space-time  formulation  was  applied  to  the  entropy-variables  formulation  to  construct 
a  fully  discrete  entropy-stable  scheme  [66].  Thereafter,  the  entropy- variables  formulation 
offered  a  foundation  for  computing  compressible  hows.  Interested  readers  are  referred  to  [35] 
for  a  detailed  review.  However,  it  needs  to  be  pointed  out  that  the  stability  of  this  entropy- 
variables  formulation  is  contingent  upon  the  convexity  of  the  entropy  function.  For  phase- 
held  problems,  the  non-convexity  of  the  entropy  function  precludes  the  possibility  of  directly 
applying  this  methodology.  To  overcome  the  challenges  posed  by  the  non-convexity  of  the 
entropy,  hrst,  the  definition  of  the  entropy  variables  is  generalized  to  the  functional  setting 
[45].  Employing  functional  entropy  variables,  we  derive  an  alternative  statement  of  the 
original  Navier-Stokes-Korteweg  equations.  The  weighted  residual  formulation  based  on  this 
alternative  statement  leads  to  a  provably  entropy-stable  semi-discrete  formulation.  Second, 
to  develop  a  stable  temporal  scheme,  we  adopt  the  methodology  based  on  special  quadrature 
rules  [24,  43,  45].  The  new  difficulty  comes  from  the  discretization  of  the  energy  time 
derivative,  since  the  isothermal  Navier-Stokes-Korteweg  equations  have  been  well  handled, 
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as  described  in  [45].  A  new  “jump”  operator  (i.e.,  discrete  time-derivative  operator)  is  devised 
for  the  total  energy.  It  is  shown  that  this  operator  represents  a  third-order  perturbation  to 
the  classical  jump  operator.  By  using  perturbed  trapezoidal  rules  repeatedly,  it  is  proven 
that  the  temporal  approximation  based  on  the  new  jump  operator  dissipates  entropy,  and  the 
requirement  of  convexity  is  released.  It  is  anticipated  that  this  new  temporal  discretization 
technology  is  applicable  to  more  general  problems. 

Non-Uniform  Rational  B-Splines  (NURBS),  in  the  setting  of  isogeometric  analysis,  are 
utilized  to  provide  a  representation  of  the  geometry  as  well  as  an  approximation  for  the 
spatial  discretization  [33].  Isogeometric  analysis  has  been  shown  to  enjoy  several  desirable 
numerical  properties:  (1)  it  retains  an  exact  representation  of  the  geometry;  (2)  it  possesses 
a  unique  fc-refinenient  technology,  which  allows  one  to  generate  higher-continuity  basis  func¬ 
tions  without  proliferation  of  degrees  of  freedom;  (3)  it  exhibits  superior  robustness  [14,  42] 
and  accuracy  [18]  properties  compared  with  traditional  finite  elements.  The  above  attributes 
make  isogeometric  analysis  a  particularly  effective  approach  in  the  approximation  of  phase- 
held  problems  [22,  25].  The  NURBS-based  initial  instantiation  of  isogeometric  analysis  has 
been  widely  used  in  both  design  and  analysis  [12].  Recent  advances  of  isogeometric  analy¬ 
sis  include  T-splines  [5,  63],  divergence-conforming  B-splines  [19,  20,  21],  and  isogeometric 
collocation  methods  [61].  In  particular,  isogeometric  collocation  methods  are  shown  to  be 
an  efficient  alternative  to  isogeometric  Galerkin  methods,  which  offers  a  potentially  powerful 
alternative  for  phase-held  simulations  [26,  59]. 

1.4  Structure  and  content  of  the  paper 

The  body  of  this  work  is  organized  as  follows.  In  Section  2,  based  on  the  Gurtin  mi¬ 
croforce  theory,  a  unihed  thermomechanical  modeling  framework  is  derived.  The  Navier- 
Stokes-Korteweg  equations,  including  the  interstitial  working  hux  term,  are  derived  within 
this  framework  by  choosing  an  appropriate  Helmholtz  free  energy  functional.  The  ther¬ 
modynamic  properties  of  this  model  are  discussed.  In  Section  3,  a  provably  entropy-stable, 
second-order  time  accurate  numerical  scheme  is  designed  and  analyzed  for  the  Navier-Stokes- 
Korteweg  equations.  Our  approach  is  based  fundamentally  on  the  concept  of  functional  en¬ 
tropy  variables  and  utilizes  newly  developed  concepts  for  the  derivation  of  time  stepping  al¬ 
gorithms  appropriate  for  problems  involving  non-convex  potentials.  In  Section  4,  benchmark 
problems  are  studied  to  verify  the  theoretical  estimates.  In  Section  5,  a  suite  of  application 
examples,  including  evaporation,  condensation,  interface  motion  under  a  temperature  gradi¬ 
ent,  and  boiling  hows,  are  numerically  investigated  using  the  model  and  algorithm  developed. 
Conclusions  and  future  research  directions  are  presented  in  Section  6. 
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2  The  Navier-Stokes-Korteweg  equations 


2.1  Balance  laws 

The  continuum  theory  setting  is  the  Euclidean  space  M3,  with  fixed  orthonormal  basis  vectors 
e*,  i  —  1,2,  3.  The  body  under  consideration  occupies  a  region  B  C  M3,  which  is  referred 
to  as  the  reference  configuration.  The  material  point  in  B  is  labeled  by  X  =  (Xi,  X2,  X^)T  ■ 
The  motion  that  the  continuum  body  undergoes  is  denoted  as  X  :  B  x  [0,  00)  — >  M3.  The 
image  of  B  by  X  at  time  t  is  denoted  as  Bt,  which  is  referred  to  as  the  current  configuration. 
The  spatial  position  of  material  points  X  at  time  t  is  given  by 


x  =  X{X,t), 


where  x  =  (xi,X2,Xs)T .  We  postulate  that  the  map  X  is  differentiable,  one-to-one,  and 
orientation  preserving  for  each  time  t  >  0.  Consider  an  arbitrary  open  set  12  of  B ;  its  image 
at  time  t  is  denoted  as  =  X  (12 ,t).  The  boundary  <912t  is  oriented  with  a  unit  outward 
normal  vector  n(x).  We  assume  that  there  exists  a  density  field  p(x,  t)  and  a  velocity  field 
u  (x,  t)  in  the  current  configuration.  The  spatial  velocity  field  is  defined  as 

d 

u  (x,  t)  —  —  X(X,  t)  where  x  =  T(X,2). 

In  the  following,  we  understand  D / Dt  to  be  the  material  time  derivative,  i.e., 


D 

Dt 


d_ 

Ft 


(■)  +  u  -  V  (■) , 


where  V  is  the  spatial  gradient  operator.  We  have  the  following  balance  laws  that  govern 
the  behavior  of  the  body. 

•  Conservation  of  Mass 


ftl^W^o.  (1) 

•  Balance  of  Linear  Momentum 

4:  [  p(x,t)u(x,t)dVx  =  f  t(x,t)dAx+  f  pb(x,t)dVx.  (2) 

Jdt  JdClt  Jilt 

Here  the  traction  field  is  given  by  t(x,  t)  =  T(x,  t)n(x),  where  T(x,  t)  is  the  Cauchy  stress 
tensor;  b(x,  t)  is  the  external  body  force  per  unit  mass. 
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•  Balance  of  Angular  Momentum 


d_ 

dt 


p(x,  t)u(x,  t)dVx 


x  t(x,  t)dAx  + 


x  b(x,  t)dVx. 


Idnt 


Ittt 


(3) 


In  our  work  the  central  modeling  subject  is  the  liquid-vapor  phase  transition.  The  phase-field 
order  parameter  for  the  change  of  the  state  of  matter  is  chosen  as  the  density  p.  Following 
the  ideas  of  Gurtin  [27],  we  assume  that  there  exists  a  set  of  forces  that  accounts  for  the 
kinematics  of  phase  transitions.  These  forces  are  called  microforces  primarily  because  they 
are  involved  with  the  local  transformation  of  the  material,  rather  than  the  macroscopic 
movements.  We  assume  that  the  kinematics  of  p  are  associated  with  the  following  forces: 

£,  the  microstress, 

X,  the  surface  microforce, 

<p,  the  internal  microforce, 

/,  the  external  microforce. 

This  set  of  microforces  is  balanced  as  stated  in  the  following  equation. 

•  Balance  of  Microforce  Associated  with  Density  Phase  Transition 


[  x(x,t)eL4x  +  f  (pdA 4+  f  ldVx  —  0.  (4) 

J  d£lt  ^  ^ 

The  surface  microforce  is  given  by  y(x,  t)  =  £(x,  t)  ■  n(x). 

Remark  1.  The  notion  of  microforce  was  initially  introduced  to  generalized  the  Cahn- 
Hilliard  equation  [27].  For  a  comprehensive  review,  interested  readers  are  referred  to  [29]. 

•  Balance  of  Energy 


jf  J  P(x,  t)E(x,  t)dVx  =  J  ^T(x,  t) u(x,  t)  +  ^P(x>  t)€(x,  t)  -  •  n(x)dAx 

/'  D 

+  b(x,t) -u(x,t) +  Z(x,t)— p(x,t) +  p(x,t)r(x,t)dyx.  (5) 

In  equation  (5),  the  following  notations  are  introduced: 

E(x,t)  =  t(x,  t)  +  ^|u(x,  t)|2,  the  total  energy  density  per  unit  mass, 

t(x,  t),  the  internal  energy  density  per  unit  mass, 

q(x,  t),  the  heat  flux, 

r(x,  t),  the  heat  source  per  unit  mass. 
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Besides  the  traditional  working  terms  of  the  macroscopic  forces  and  the  macroscopic  sources, 
there  are  non-classical  terms  contributing  to  the  change  of  the  total  energy.  These  terms  are 
the  power  expenditures  of  the  microstress  £  and  the  external  microforce  /: 


D 

Dt 


p(x,t)£(x,t)  •  n(x)oL4x, 


'fit 


The  internal  microforce  does  not  contribute  to  the  energy  change.  See  [27]  for  a  conceptual 
explanation. 

•  Second  Law  of  Thermodynamics 


d  r 

V(x,t)dVx  :=—  /  p(x,t)s(x,t)dVx  + 


'fit 


dt 


'fit 


'd  fit 


q CM)  •  n(x) 
d(x,t) 


dAy 


P(x,t)r(x,f) 

P(x,t) 


filC  >  0. 


(6) 


Here  T>(x,t)  denotes  the  total  dissipation,  s(x,t)  denotes  the  entropy  density,  and  9{x,t )  is 
the  absolute  temperature.  This  inequality  is  the  second  law  of  thermodynamics,  or  Clausius- 
Duhem  inequality. 

Applying  the  divergence  and  the  Reynolds’  transport  theorems,  we  can  obtain  the  gov¬ 
erning  equations  and  the  Clausius-Duhem  inequality  in  local  forms  (omitting  the  arguments 
x  and  t  for  simplicity)  as 


TTt+pV' u  =  0’ 

(7) 

plV  =  VT  +  pb' 

(8) 

T  =  Tt, 

(9) 

V  •  £  +  (f  +  l  =  0, 

(10) 

DE  (  Dp  \  Dp 

pDt=v'[Tu+Dd  q)+pb'u+Tt+pr’ 

(11) 

^  Ds  „  /q\  pr 

v-=pm+vil)-J>-a- 

(12) 

In  addition  to  the  governing  equations  (T)— (12) ,  there  is  a  balance  equation  for  the 
internal  energy  i.  The  energy  balance  equation  (11)  can  be  expanded  as 


Dl  Du  Dp  /  Dp 

pm+pu' m  =v  t.u  +  t:Vu  +  v.«— +  «  v(— 

_  ,  .Dp 

-  V  ■  q  +  pb  •  u  +  /  —  +  pr. 


(13) 


The  linear  momentum  balance  equations  and  the  microforce  balance  equation  can  be  utilized 
to  give  the  following  relations. 


Du 

pu  ■  -jy-  =  V  •  T  •  u  +  pb  •  u, 
Y7  eDp  <iDp  Dp 

^  Dt  Dt  Dt 


(14) 

(15) 


Substituting  (14)-(15)  into  (13),  we  may  obtain  a  balance  equation  for  the  internal  energy 
as  follows. 


Dt  Dp  /  Dp' 

pWt  =  T:Vu-v-a+i-v{m^-v'q  +  pr- 


(16) 


This  equation  will  be  used  as  a  starting  point  for  the  derivation  of  constitutive  relations. 


2.2  Coleman-Noll  type  analysis  and  constitutive  relations 

To  close  the  model,  we  still  need  to  provide  the  constitutive  relations  for  the  Cauchy  stress, 
the  internal  energy  density,  the  entropy  density,  the  heat  flux,  and  the  microforces.  In  this 
section,  we  derive  the  explicit  form  of  the  constitutive  relations  in  terms  of  a  thermodynamic 
potential.  In  this  derivation,  the  Coleman-Noll  argument  [11]  is  applied  so  that  the  resulting 
constitutive  relations  will  be  thermodynamically  consistent. 


2.2.1  Free  energy  imbalance 

The  Helmholtz  free  energy  density  per  unit  mass  T  (x,  t)  is  defined  by 

T(x, t)  :=  t(x, t)  —  #(x, t)s(x, t). 

Taking  material  time  derivatives  at  both  sides,  we  get  the  relation 

Dl  Ds  D T  DO 

- 6 —  = - b  s — . 

Dt  Dt  Dt  Dt 


(17) 
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Substituting  the  internal  energy  balance  equation  (16)  and  the  second  law  of  thermodynamics 
(12)  into  the  above  relation,  we  can  get  an  inequality 


DV  DO  Dp  f  Dp\  q-W 

dwrr  +  ps-r  <  T  :  Vu  -  <p-L  +  $.  •  V  — 


Dt  Dt 


Dt 


Dt 


e 


Moving  the  term  psDO / Dt  to  the  right  hand  side,  we  can  get  a  constraint  inequality  for  T 
as 


DT  m  „  Dp  ^  f  Dp\  q-W  DO 
P^r  <  T  :  Vu  -  <p-L  +  £  •  V  -M - — - ps- 


Dt 


Dt 


Dt 


0 


Dt 


(18) 


The  inequality  (18)  is  referred  to  as  the  free  energy  imbalance.  It  plays  an  analogous  role 
to  (12)  in  placing  restrictions  on  constitutive  relations.  In  fact,  for  pure  mechanical  pro¬ 
cesses  when  thermal  effects  are  negligible,  the  Helmholtz  free  energy  is  the  thermodynamic 
potential  that  characterizes  the  dissipation  behavior  of  the  isothermal  system  [62],  For  an 
isobaric  isothermal  process,  however,  the  Gibbs  free  energy  should  be  chosen  as  a  proper 
thermodynamic  potential  [46,  62].  In  this  work,  an  isobaric  process  is  not  assumed.  Hence, 
the  Helmholtz  free  energy  is  a  valid  thermodynamic  potential.  Before  proceeding  further, 
we  split  the  Cauchy  stress  T  and  the  velocity  gradient  Vu  into  deviatoric  and  hydrostatic 
parts. 


1.  The  Cauchy  stress  T  can  be  split  into  deviatoric  and  hydrostatic  parts, 


T 


rpd  _|_  rji/l 


(19) 


where 

Td  =  T  -  -  (trT)  I, 

3 

Th  =  ^  (trT)  I. 

Here  I  is  the  identity  tensor,  and  tr(-)  is  the  trace  operator. 
2.  The  velocity  gradient  can  be  split  into  three  parts, 

Vu  =  Ld  +  Lh  +  W, 


(20) 

(21) 


(22) 


wherein 


Ld  =  -  (Vu  +  VuT)  -  -V  •  ul, 

2  3 


(23) 
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(24) 

(25) 


Lh  =  -V  •  ul. 

3 

W  =  ^  (Vu  -  VuT)  . 

In  this  split.  1/  and  Lh  are  the  deviatoric  and  hydrostatic  (i.e.,  dilatational)  parts  of 
the  rate  of  strain  tensor  L;  W  is  the  spin  tensor  (i.e.,  vorticity  tensor). 

Consequently,  it  is  straightforward  to  make  the  following  observations. 

1.  According  to  the  mass  conservation  equation  (7),  we  have 

V-u  (26) 

P 

2.  The  gradient  of  material  time  derivative  Dp/ Dt  can  be  expanded  as 

v  (§? ) =  <Vp> + VuTVp  =  §i (Vp) + L''Vp + wTVp  ~  ^r~Vp-  (27) 


3.  Making  use  of  the  property  of  deviatoric  tensors,  the  inner  product  of  T  and  Vu  can 
be  written  alternatively  as 

T  :  Vu  =  Td  :  Ld  +  Th  :  Lh  =  Td  :  Lcl  +  -(trT)V  ■  u.  (28) 

3 


Making  use  of  the  above  observations,  the  free  energy  imbalance  (18)  relation  can  be  rewrit¬ 
ten  as 


> - 

Dt 


<T< 


'  :Ld  — 

q  -we 
9 


trT  Dp  Dp  D  .  , 

yrK“^+«w(Vp)+Vp'L£+VpW£ 


ps 


DO 

Dt' 


1 

3 ~p 


X7p-£ 


Dp 

Dt 

(29) 


2.2.2  Coleman-Noll  type  analysis 

Invoking  Truesdell’s  principle  of  equipresence  [70],  we  assume  that  T,  s,  T,  £,  q,  i,  and  tp 
are  functions  depending  on  p,  Vp,  Dp/Dt ,  9,  V0,  and  LJ.  Specifically,  the  Helmholtz  free 
energy  density  T  can  be  written  as 

,e,v«.iA 


T  =  T 


„  Dp 

pyp’m 
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The  material  time  derivative  of  T  and  the  chain  rule  lead  to 


DT_«9TDp  |  <9T  D(Vp)  ,  <9T  D2p  D9 

~Dt  ~~dp~Dt  +  d  {Vp)  '  Dt  +  c>  (Dp/Dt)  DD  +  +  0(V0) 

«9T  DLd 

+  '  ~W' 


d{vo) 

Dt 


(30) 


Now  substituting  (30)  into  the  free  energy  imbalance  (29)  and  making  use  of  the  relations 
(22)-(28),  we  can  get 


P 


DT 

Dt 


<9T  Dp  <9T  D  {Vp) 


—P\  + 


chi/ 


dp  Dt  '  d(Vp) 
<9T  DLfi 


Dt 


D2p  <9T  D0 

+  9  (Dp/Dt)  Dt2  +  Dt  +  c>  (V0) 


<9T  D  (V0) 


Dt 


<9IH  '  Dt 


<  Td  :  Ld 


trT  Dp  Dp  D  , 

yr«^iv+^iM(V'>)+v'’'L« 


„  ,,,,  1  Dp  qV(  £W 

+  Vp.w«--V^«— -  — 

Grouping  terms  together,  the  above  inequality  is  reorganized  as 

chi/ 


<9T  trT 

chi/ 


*  +  £v'-«)^+W(vw 


D0 


<9T 


p¥¥ +  psJ  Dt+pd(W) 


—  W  :  (Vp  <8>  £)  <  0. 


D  (V0)  q  ■  V0  <9T 

+  -L^—  +  p 


DL° 


(31) 


D  <9T  D2p 

Dt  ^  +  P~d{D^fDt)~DD 


Dt 
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3Ld  ‘  Dt 


-  1/  :  (Td  +  Vp  (8)  |) 
(32) 


Here  we  provide  an  analysis  of  (32)  by  invoking  the  arguments  made  by  Coleman  and  Noll 
[11].  We  notice  that  (32)  is  linear  in  D(Vp)/Dt,  D2p/Dt2,  DQ/Dt,  D(V0)/Dt,  DLC,/Dt, 
and  W.  Through  appropriate  choices  of  external  forces  and  external  sources,  we  may  have 
arbitrary  levels  of  the  material  rates  of  the  state  variables  and  the  spin  tensor  in  (32)  at  a 
particular  time.  Hence,  in  order  to  satisfy  (32),  we  must  have 


&  <9T 

€~pa(vPy 

(33) 

(34) 

d  {Dp/Dt)  ~  °’ 

chp 

(35) 

s  = 

chp  ^ 

(36) 

d{V9)  _°’ 
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dV 

d\? 


chi' 

dUl 


T 


(37) 


Vpg)£=£®Vp. 


(38) 


Due  to  the  symmetry  of  L(/,  d^/DL'1  has  to  be  symmetric.  This  fact  together  with  (37) 
implies 


<9T 

d\? 


0. 


(39) 


According  to  (34),  (36),  and  (39),  the  Helmholtz  free  energy  density  T  has  to  be  independent 
of  Dp/Dt ,  VO ,  and  L(/.  It  can  be  written  as 


T  =  T(p,Vp,0). 


(40) 


The  relation  (32)  is  reduced  to 


+ 


trT 


+  <p  + 


Dp  q  ■  VO 
~Dt  +  0 


Ld  .  (Td  +  Vp  <g>  £)  <  0. 


(41) 


Based  on  the  above  inequality,  the  following  choices  are  made.  While  these  relations  are  not 
the  most  general  possible  forms  that  can  satisfy  (41),  we  will  show  that  they  do  suffice  to 
recover  the  most  common  physical  models. 


trT  <9T  1  £>p 

*  =  “  "sW  *  ct  ’ 

q  =  —kVO, 

TJ  =  2pLJ  -  i  (Vp  +  (  ®  Vp)  +  tVp  ■  «I. 


(42) 

(43) 

(44) 


where  B,  k,  and  p.  are  scalar  valued  functions. 

Remark  2.  The  relation  (38)  can  be  viewed  as  a  result  of  objectivity.  In  conjunction  with 
(33),  it  can  be  rewritten  as 


Vp®  p 


dV 

0(Vp) 


P0(Vp) 


(8)  Vp. 


(45) 


It  poses  a  constraint  on  how  the  non-local  gradient  term  Vp  can  enter  into  the  free  energy 
density  function  T.  It  has  been  shown  in  [9]  that  terms  in  the  form  of  t)  •  Vp  or  €  :  V2p  , 
with  h  being  a  constant  vector  and  €  being  a  constant  second- order  tensor,  cannot  enter  into 
the  free  energy  density  function  T.  This  assertion  can  be  easily  justified  using  the  constraint 
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relation  (45).  We  will  only  consider  the  case  with  Vp  appearing  in  the  free  energy  density 
function  4/  as  |Vp|2.  This  is  in  fact  the  case  considered  by  van  der  Waals  in  his  seminal 
work  [58],  This  special  choice  for  the  Helmholtz  free  energy  density  function  guarantees  the 
satisfaction  of  the  relation  (38) . 

Remark  3.  It  should  be  pointed  out  that  (42) -(44)  are  specific  choices  made  to  guarantee  the 
inequality  (41).  It  is  feasible  to  propose  more  general  constitutive  relations  for  <p,  q,  and  Tf/ 
which  satisfy  the  inequality  (41).  However,  discussion  of  such  a  general  constitutive  theory 
is  beyond  the  scope  of  this  work. 

2.2.3  Constitutive  relations 

Based  on  the  relations  (42)- (44),  we  can  obtain  the  constitutive  relations  expressed  in  terms 
of  T. 

Microstress 

The  relation  (33)  gives  the  constitutive  relation  for  the  microstress, 

t  _ 

i~pdWp) 


(46) 


Cauchy  stress 

From  the  relations  (33),  (42),  (43),  and  the  microforce  balance  equation  (10),  we  have 
the  constitutive  relation  for  the  hydrostatic  part  of  the  Cauchy  stress  as 


T 


h 


/  <9T  \  2 ST 


!  d T 
3  P<9(Vp) 


■  Vp  +  pi  —  Bp 


I. 


(47) 


Replacing  the  microstress  £  by  the  relation  (46),  the  deviatoric  part  of  the  Cauchy 
stress  is  given  by  (44). 


T(/  =2pLd 


£(vp< 


<9T 


+ 


<9T 


<9(Vp)  d  (Vp) 


Vp 


d(V,i) 


i 


(48) 


Combining  the  two  parts,  the  Cauchy  stress  T  reads  as 


rp  _ rpd  _|_  rpAi 


—2fiLd  -  P-  [  Vp 


ST 


+ 


ST 


d(Vp)  d  (Vp) 


Vp 


pV-  p 


,<9T 


Dp 


d  (Vp)  J  dp 


+  pi  —  Bp—T-  I. 


Dt 


(49) 
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Heat  flux 


The  constitutive  relation  for  the  heat  flux,  given  by  (43),  is  the  Fourier’s  law  [48]: 

q  =  —nV6,  (50) 

wherein  k  is  the  thermal  conductivity. 

Entropy 

The  relation  (35)  defines  the  entropy  density  s  as 

s  =  ~W 

This  definition  coincides  with  the  classical  thermodynamic  definition  [48] 
quently,  the  internal  energy  density  i  is  given  by 

i  =  V  +  Os  =  T  -  9— . 

o0 

2.3  Dissipation  inequalities 

In  this  section,  the  choices  (42)-(44)  are  validated  by  analyzing  the  dissipation  of  the  model. 
It  will  be  clear  how  different  terms  enter  into  the  dissipative  mechanisms  in  isolated  and 
isothermal  processes. 

Lemma  1.  Given  the  constitutive  relations  (46) -(52),  the  dissipation  V  defined  in  (12)  takes 
the  form 

vGfplLr  +  lB^y  +  ±Km>.  (53) 

Proof.  We  start  by  considering  the  internal  energy  balance  equation 

4  =  T:Vu-^+«'v©-v'q+'”''  (54) 

It  is  known  from  (28)  that 

T  :  Vu  =  Td  :  L(Z  +  -  (trT)  V  •  u.  (55) 

3 


(51) 
Conse- 

(52) 
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Making  use  of  the  constitutive  relation  (48),  we  have 


T<! .  I/!  =Ld  .  2jXLd  -  J  |  Vp 
=2fi\L‘\ 2-|(Vp 


<94/  <94/ 

+ 


(9  (Vp)  '  <9  (' Vp)  ~  VP)  ''  L  +  3Vp'«9(Vp) 


<94/ 


d,  ,  P  V7  „  J  . 


<94/  <94^ 

+ 


<9  (Vp)  <9  ( Vp) 
According  to  the  constitutive  relation  (47),  we  have 


Vp  :  L" 


(56) 


-  (trT )  V  •  u  =pV  ■  (  p 


\  29vi/  i  9$ 

^)Jv-u-p^v'u-3p5(v^rVpV'11 


+  p/V  •  u  —  Rp^^V  •  u 

/  94  \  Dp  <94<  Dp  1  94  Dp 

'  \Pd  (Vp) )  ~Dt  +  +  3  9  (Vp)  ‘ 

+Sp2(v'u)2' 


(57) 


Recalling  from  the  relation  (27),  we  have 


<-viW>  =< 


=p- 


4)  (Vp)  +  L"vp  +  wTvp  -  Vp 


<94/ 


4)  (Vp)  +  LWp  +  wrvp  -  444H  Vp 


d  (Vp) 

The  microforce  balance  equation  implies  <p  =  —V  ■  £  —  l.  Consequently,  we  have 


(58) 


DP  „  ,-Dp  ,-Dp 
y  Dt  *  Dt  dD 


(59) 


Now  substituting  (55)- (59)  into  (54),  and  using  (38)  repeatedly,  we  obtain 


p^=2p|L"|2-^  (Vp. 


<9T  <9T 

+ 


<9  (Vp)  <9  (Vp) 


Vp);L,‘-vipafy)^ 


P 


P 


d^Dp  1  94  Dp  Dp  2  a2  ,  ^  &Dp  ,  ,Dp 
tyDi  3WP)  PDi  ~lDi  +  Bp  (V'U)  +V'l'K+,'« 


(9T 

a(vp) 


4)  (vp)  +  L"vp  +  wTvp  -  Vp 


V  ■  q  +  pr 


=  2p\L'l\~  +  Bp2  (V  •  u)2  +  p9^-1^-  +  Py^y  •  ~~  (Vp)  -  V  ■  q  +  pr 


dp  Dt  <9  (Vp)  Dt 
Z9T  <9TL>9\ 


^2p|Lr  +  V(V.ur  +  p(— ---j-V.q  +  pr 
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—  V  ■  q  +  pr. 


-  2£|Lf  +  Bp2  (V  ■  u)2  +  P  (j|  - 


Moving  all  time  derivative  terms  to  the  left  hand  side  yields 

l2  +  iV(V.u)2-iv.q+ipr. 


By  definition,  we  have 


D:=p 

_1 

_1 


Ds 


VV7  +  V-  ^ 


Dt 


q 


pr 


e 


2/f|I/f  +  -Hp2  (V  •  u)2  - 

u 


2p\Ld\2  + -Bp2  (V  -u)2  + 

u 


q  •  V0 

d2 


1 

d2 


«|Vd|2, 


which  completes  the  proof.  □ 

The  dissipation  formulation  (53)  suggests  that  the  model  will  guarantee  the  second  law 
of  thermodynamics  if  the  material  moduli  are  non-negative,  which  is  summarized  in  the 
following  theorem. 

Theorem  1.  If  ft,  k,  and  B  are  non-negative,  the  system  of  balance  equations  (7)-(ll) 
satisfies  the  second  law  of  thermodynamics  in  the  following  sense. 

V  =  t2p|L-'|2  +  tflp2  (V  ■  u)2  +  iK|V«|2  >  0, 


The  proof  of  this  theorem  follows  straightforwardly  from  Lemma  1.  The  significance  of 
this  theorem  is  that  the  modeler  only  needs  to  design  an  explicit  formulation  for  the  thermo¬ 
dynamic  potential.  Once  it  is  given,  the  model  is  closed  with  non- negative  dissipation.  Under 
isothermal  conditions,  the  entropy  dissipation  relation  will  degenerate  into  an  inequality  for 
the  summation  of  the  Helmholtz  free  energy  and  the  kinetic  energy. 

Lemma  2.  Under  the  isothermal  condition,  if  u  =  0  and  £  •  n  =  0  on  the  boundary  dQt,  the 
following  relation  holds. 


d_ 

dt 


dUx 


•  u  +  l 


Dp 

Dt 


dVx. 


(60) 


Proof.  Since  6  is  constant,  according  to  T  =  l  —  6s,  one  has 


D T  _  Dl  yDs 
Dt  Dt  Dt 
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and 


1  i 

T  H — lul2  =  l  -\ — |u 1 2  —  9s. 

2  2 

Multiplying  the  the  above  equation  by  p  and  integrating  over  Tlt  result  in 

it  L. p  (® + ^ |u|2)  214 =i  L. p  (* + ^ |u|2  ■ 0s ) 2,4 

=L(^+^)-"dA*+L(pb-u+i%-9v)dv'- 

The  boundary  integral  terms  are  canceled  due  to  the  boundary  conditions,  and  hence 

i  L p  (® + ^|u|2)  ^ = L  (pb ■ u + 1%  - ev) dv" 

which  completes  the  proof  of  the  lemma.  □ 

Based  on  Lemma  2,  we  may  obtain  the  following  stability  theorem  for  isothermal  pro¬ 
cesses. 

Theorem  2.  If  (1)  the  system  undergoes  an  isothermal  process,  (2)  u  =  0  and  £  •  n  =  0  on 
the  boundary  d Ot ,  (3)  the  forces  b  =  0  and  l  —  0  in  Qt,  and  (4)  the  material  moduli  p  and 
B  are  non-negative,  the  stability  of  the  system  is  given  by  the  following  dissipation  relation. 

jt  £  P  (*  +  dV*  =  -  fQ  (2^ILT  +  BP2  (V  '  u)2)  dV *  <  °- 

Remark  4.  According  to  the  constitutive  relation  (46)  ,  the  boundary  condition  £  ■  n  =  0  is 
equivalent  to  Vp  •  n  =  0.  The  general  contact-angle  boundary  condition  is 

~m\ ' 11  =  cosW’  (61) 

wherein 


IIWII  =  vWp.vp, 

mid  cf)  is  the  contact  angle  of  the  diffuse-interface  against  the  wall  boundary  measured  in  the 
vapor  phase  (see  Figure  1).  Hence,  Vp-n  =  0  gives  the  ninety-degree  contact  angle  boundary 
condition. 
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Figure  1:  Illustration  of  the  contact  angle  boundary  condition  (61).  The  red  arrow  points  in 
the  direction  — Vp/||Vp||. 

2.4  The  van  der  Waals  fluid  model 

In  the  preceding  section,  a  general  continuum  mechanics  modeling  framework  has  been  estab¬ 
lished,  with  the  objective  of  taking  non-local  effects  into  account.  Theorems  1  and  2  reveal 
that  the  model  is  thermodynamically  consistent  if  the  material  moduli  are  non-negative. 
Thus,  the  modeling  work  is  principally  reduced  to  a  proper  design  of  the  thermodynamic 
potential.  This  design  procedure  is  primarily  based  on  the  consideration  of  thermodynamics. 
Our  discussion  will  focus  on  the  van  der  Waals  fluid.  The  full  thermomechanical  theory  of 
the  van  der  Waals  fluid,  initially  derived  by  Dunn  and  Serrin  [17],  will  be  recovered.  We  will 
discuss  preliminary  thermodynamic  properties  of  the  system. 

2.4.1  Governing  equations 

Van  der  Waals’  theory  [58]  is  considered  well-suited  for  describing  liquid-vapor  phase  tran¬ 
sitions.  In  thermodynamics,  the  Helmholtz  free  energy  density  for  the  van  der  Waals  fluid, 
T.  is  given  by 

Up,0,Vp)  =  4V(p,0)  +  Vvp|2,  (62) 

=  -ap  +  RO log  (^  P_  I  -  CVB log  (T-J  +  CJ,  (63) 

where  a,  b  are  associated  with  fluid  properties  whose  meanings  will  be  revealed  in  the  coming 
discussion;  9ref  >  0  is  the  reference  temperature  for  the  model;  R  is  the  specific  gas  constant; 
Cv  is  the  specific  heat  capacity  at  constant  volume  for  the  van  der  Waals  fluid;  A  is  the 
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capillarity  coefficient.  In  this  work,  we  assume  A  is  a  constant.  With  the  Helmholtz  free 
energy  given,  the  constitutive  relations  can  be  readily  obtained.  According  to  (46),  the 
microstress  for  the  van  der  Waals  fluid  is 


chi' 

*  =  =  AVp. 

d  (Vp) 

Following  (49),  the  Cauchy  stress  can  be  written  explicitly  as 


m  T  .  p  f  (9'h  9'h 

=  ^  ~~  2  (Vp  s  Wjypj +  &Wp) 


Vp 


pv  •  i  p  _ ,  -p2-^  +  p/  +  Hp2v-u  i 


<94'  A  2  <94/ 

J<9  (Vp) ;  _p  ~dp 

=2pLd  -  AVp  (8)  Vp  +  (^ApAp-p2^^  +  ^|Vp|2  +  p/  +  Hp2V-u  )  I. 
For  convenience,  the  Cauchy  stress  can  be  split  into  three  parts: 


T  =  t  +  q  -  pi, 


(64) 


wherein 


t  —2pLd  +  Bp2V  ■  ul, 


A, 


=  —  AVp  <8>  Vp  +  |^ApAp  +  V |  Vp|2  +  pi  )  I, 

2^^loc  £,7  P^  2 

P  =P  -w —  =  /to- - op  . 


<9p 


b  —  p 


(65) 

(66) 

(67) 


Here,  r  represents  the  viscous  stress,  £  represents  the  Korteweg  stress,  and  p  stands  for  the 
thermodynamic  pressure.  In  the  subsequent  discussion,  we  assume 


in  which  p  and  A  are  the  first  and  second  viscosity  coefficients  [56].  With  these  choices, 

<?  =  -  AVp  <8)  Vp  +  f  ApAp  +  Vp|2^  I, 
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which  is  identical  to  the  stress  derived  by  and  named  after  Korteweg  [39],  and  r  is  the 
viscous  stress  for  Newtonian  fluids: 


r  =/j  (Vu  +  VuT)  +  AV  •  ul. 

The  entropy  density  s  and  the  internal  energy  density  i  are 

s  =  -  i?  log  +  C„  log  (£-)  , 

i  =  —  ap  +  Cv9  +  —  |  V  p  | 2 . 

Zp 

We  recall  that  in  the  derivation  of  the  constitutive  relation  for  T/;  in  (47),  the  microforce 
balance  equation  is  used.  Similarly  to  the  angular  momentum  balance  equation,  the  micro¬ 
force  balance  equation  (10)  is  satisfied  by  the  constitutive  relations  and  decoupled  from  the 
system.  Let  us  denote  the  power  expenditure  of  the  microstress  as 

II  =  £Dp/Dt  =  ApV  ■  uVp. 

The  governing  equations  for  the  van  der  Waals  fluid  are 

^  +  V  •  (pu)  =  0,  (68) 

+  V  •  (pu  <g>  u)  +  Vp  —  V  •  r  —  V  •  <?  =  pb,  (69) 

V  '  ((PE  +  p)  u  -  (r  +  q)  u)  +  V  ■  q  +  V  •  n  =  pb  •  u  +  pr.  (70) 

This  system  of  equations  is  known  as  the  Navier-Stokes-Korteweg  equations  [17].  According 
to  Lemma  1,  the  dissipation  for  the  system  is 

V  =  ^IL"!2  +  l-Bp2  (V  ■  u)'2  +  iK|V0|2 

=  ^2p|L"|2  +  l  (*  +  V  •  u)2  +  L/c|  V9|2.  (71) 

To  ensure  the  second  law  of  thermodynamics,  it  is  sufficient  to  require  that 

—  2  _ 

p  >  0,  A  +  -p  >  0,  k  >  0. 

o 

The  term  II  was  initially  introduced  by  Dunn  and  Serrin  to  enforce  the  thermodynamic 
consistency  and  was  referred  to  as  the  “interstitial  working  flux”  [17].  In  our  framework, 
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V  •  II  appears  naturally  as  the  power  expenditure  of  the  microstress  The  previously 
mysterious  term  finds  a  rational  mechanics  explanation  in  the  microforce  theory  [29]. 

Remark  5.  The  viscous  stress  r  can  be  reorganized  as 


t  =  2//.L  +  AV  •  ul 


=  2pLd  +  2pLh  +  AV  •  ul 


=  2pLd  + 


Y  2  _ 


V  •  ul. 


Making  use  of  the  split  o/Vu  given  in  (22),  one  has 

t  :  Vu  =  (2 ft Ld  +  ^A  +  ^  V  •  ul ^  :  (L'z  +  Lh  +  W) 

—  2p|L(/|2  +  ^A  +  -p'j  (V  •  u y  . 

Hence,  the  dissipation  V  can  be  written  equivalently  as 

»  =  ir;Vu+t;K|Ve|2.  (72) 

Remark  6.  Recalling  that  the  capillarity  coefficient  A  is  assumed  to  be  constant,  the  term 
V  •  can  be  written  in  the  following  non-conservative  form. 


V  •  c  =  ApV  (Ap) .  (73) 

Remark  7.  Choosing  the  Helmholtz  free  energy  density  function  as 

*  =  RB  log  (p)  -  cve  log  ( -M  +  cve, 

V  a  ref  J 

the  compressible  Navier-Stokes  equations  can  be  recovered. 


2.4.2  Thermodynamic  properties 

We  start  the  discussion  on  the  thermodynamic  properties  by  defining  the  critical  point.  The 
critical  point  ( pCrit,&crit )  is  defined  to  be  the  values  of  density  and  temperature  that  satisfy 


dp 

dp 


(P 


criti 


0, 


d2p 
dp 2 


0. 
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Figure  2:  Illustration  of  the  van  der  Waals  pressure  p  given  by  (67)  at  different  temperatures. 
The  colored  squares  delimit  the  elliptic  regions.  The  critical  point  is  marked  by  a  black  circle. 


Simple  calculations  show  that  the  solutions  of  the  above  equations  are 


Pcrit 


b 

3’ 


6 


crit 


8  ab 
27R' 


and  the  critical  pressure  is  pcrit  :=  p(pcrit,  6 crit)  =  ab2 /27.  In  Figure  2,  the  van  der  Waals  pres¬ 
sure  function  is  plotted  as  a  function  of  density  by  fixing  the  temperature.  It  can  be  observed 
from  the  figure  that  the  pressure  function  is  not  monotone  when  the  temperature  is  below  the 
critical  temperature,  and  there  is  a  region  in  which  the  pressure  decreases  with  the  increase 
of  the  density.  This  region  is  commonly  referred  to  as  the  elliptic  region,  since  the  system 
of  balance  equations  is  of  the  first-order  elliptic  type  in  the  vanishing  viscosity-capillarity 
limit  within  this  region.  The  approximation  property  of  the  van  der  Waals  equation  of 
state  is  demonstrated  by  comparing  with  the  data  for  real  fluids.  In  Figure  3,  the  van  der 
Waals  equation  of  state  is  plotted  as  a  function  of  density  at  temperature  6  =  0.85 9crit  (blue 
solid  curve).  The  thermodynamic  data  for  water,  carbon  dioxide,  methane,  and  propane  are 
downloaded  from  the  NIST  database  [51]  and  plotted  in  the  same  figure  in  dimensionless 
form.  As  can  be  seen,  the  van  der  Waals  model  gives  a  qualitatively  accurate  description  of 
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(a) 


(b) 


Figure  3:  Comparison  of  the  van  der  Waals  equation  of  state  with  real  fluids  at  temperature 
9  =  0.85 9crit.  The  data  for  water,  carbon  dioxide,  methane,  and  propane  are  obtained  from 
[51]  and  scaled  to  dimensionless  form.  Figure  (b)  gives  a  detailed  view  in  the  vapor  phase. 
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various  fluids  in  both  vapor  and  liquid  states.  Considering  a  binomial  expansion 


('-Q  'k1+s+f’  when|?|«1’ 

the  thermodynamic  pressure  can  be  approximated  as 

+  (74) 

when  the  density  is  small.  This  suggests  that  the  van  der  Waals  theory  can  be  viewed  as 
a  high-order  modification  of  the  perfect  gas  law.  Nowadays,  modifications  of  the  van  der 
Waals  model  are  introduced  by  adding  more  high-order  terms  to  improve  the  approximation 
property  for  specific  materials.  Examples  include  the  Beattie-Bridgeman  equation  [6]  and 
the  Benedict- Webb- Rubin  equation  [7].  Another  modification  was  made  by  Serrin  [65],  who 
introduced  a  new  equation  of  state  in  the  form 

pserrin  =  Rb-^-  -  a9~sp\ 

b  —  p 

wherein  s  <  1  and  r  >  1  are  two  parameters.  This  model  was  claimed  to  give  very  accurate 
pressures  over  a  large  range  of  temperature. 

Next,  let  us  introduce  the  local  electrochemical  potential  zyoc  as 

_<9(pT/oc) 

I^IOC  *  r\ 

op 

It  does  not  come  from  our  preceding  thermomechanical  theory.  It  is  a  purely  thermodynamic 
quantity.  With  the  local  Helmholtz  free  energy  function  given  in  (63),  the  electrochemical 
potential  can  be  written  explicitly  as 

uloc  =  -2  ap  +  R6  log  +^~  loS  (/^)  +  C*e- 

The  equilibrium  state  at  a  given  temperature  can  be  determined  by  constructing  a  com¬ 
mon  tangent  line  passing  thorough  the  free  energy  curve  pflhoc  at  two  points  {pu  pi^ ioc{pi)) 
and  (p„,  Pv^ioc(Pv))-  These  two  points  correspond  to  the  energetically  stable  liquid  and  vapor 
states  at  the  temperature  and  are  usually  referred  to  as  the  Maxwell  states.  Mathematically, 
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Figure  4:  Illustration  of  the  local  free  energy  p4boc  of  the  van  der  Waals  fluid  given  by  (63) 
at  temperature  6  =  0.8 9crit-  The  green  squares  delimit  the  elliptic  region.  The  magenta 
dash-dot-line  is  the  common  tangent  line  passing  thorough  the  Maxwell  states,  which  are 
marked  as  the  magenta  circles. 


the  common  tangent  line  requires  that 


<9(pTioc)  <9(p4qoc) 

\Pv)  -  — x - (Pz), 


dp 


dp 


d(p^loc)  ,  n  ^  „  lTr  I  „  \  _  „  ^(pTzoc)  ( n\  „lTr  N 

pv  rj  ^  \pv )  pv^x  locypv)  pi  ^  \pl)  Pl^loC\Pl)‘ 


dp 


dp 


(75) 

(76) 


The  relation  (75)  implies  the  local  electrochemical  potentials  vioc  at  the  two  states  are  iden¬ 
tical.  The  relation  (76)  can  be  rewritten  as 


P 


2  d^loc  [  ^  2  d^  loc  f 

{pv)  =  Pi—^r{pih 


dp 


dp 


or,  equivalently,  p(pv)  =  p{pi)-  Therefore,  the  system  is  in  electrochemical  and  mechanical 
equilibrium  at  the  Maxwell  states.  The  Maxwell  states  together  with  the  common  tangent 
line  are  illustrated  in  Figure  4.  It  can  be  clearly  observed  from  the  figure  that  the  common 
tangent  line  lies  below  the  energy  curve,  which  implies  the  two-phase  state  is  favored  against 
the  homogeneous  mixture  state,  according  to  the  minimum  energy  principle. 

The  thermodynamic  properties  of  the  van  der  Waals  fluid  model  can  be  better  understood 
by  drawing  a  6-p  phase  diagram.  In  Figure  5,  the  elliptic  region  is  circumscribed  by  the 
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Figure  5:  Illustration  of  the  elliptic  region,  the  metastable  regions,  the  spinodal  line,  and 
the  binodal  line  for  the  van  der  Waals  fluid. 


dashed  spinodal  line  and  is  colored  in  grey.  By  connecting  the  Maxwell  states,  we  get  the 
binodal  line,  which  is  drawn  as  the  black  solid  curve  in  Figure  5.  The  regions  enclosed  by 
the  binodal  line  and  the  spindoal  line  are  the  metastable  liquid  and  vapor  regions,  which  are 
colored  in  green  and  blue  respectively.  The  metastable  states  are  physically  accessible  but 
energetically  unstable  [37].  With  enough  thermodynamic  perturbations,  the  energy  barrier 
may  be  overcome  and  the  metastable  states  may  evolve  toward  a  more  stable  two-phase 
system.  The  binodal  line  and  the  spinodal  line  meet  at  the  critical  point.  Above  the  critical 
temperature,  the  fluid  becomes  supercritical,  and  there  are  no  more  distinct  liquid  and  vapor 
states  [62], 

3  Numerical  analysis 

In  this  section,  we  focus  on  the  design  of  numerical  schemes  for  the  Navier-Stokes-Korteweg 
equations  that  preserve  the  dissipation  property  of  the  original  strong-form  problem. 
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3.1  Initial-boundary  value  problem  for  the  Navier-Stokes-Korteweg 
equations 

We  consider  a  fixed,  open,  connected,  and  bounded  domain  12  C  M3.  The  boundary  of  12 
is  denoted  as  <912  and  is  assumed  to  be  sufficiently  smooth.  The  time  interval  of  interest  is 
denoted  (0 ,  T),  with  T  >  0.  The  Navier-Stokes-Korteweg  equations  are  considered  in  the 
space-time  domain  12  x  (0,  T)  as 

^  +  V  •  (pu)  =  0,  (77) 

+  V  •  (pu  ®  u  +  pi)  -  V  •  r  -  V  •  <;  =  pb,  (78) 

d  ( oE) 

— +  V  •  (( pE  +  p)  u  -  (r  +  <j)u)  +  V  ■  q  +  V  •  II  =  pb  •  u  +  pr.  (79) 

In  this  section,  we  impose  periodic  boundary  conditions  for  all  variables.  Therefore,  the 
problem  can  be  regarded  to  be  a  periodic  flow  posed  on  a  three-dimensional  torus  T3  in 
space.  Given  p0  :  12  — »  (0,  b),  u0  :  12  — >  Md,  and  6*0  :  12  — >  M  as  the  initial  density,  velocity, 
and  temperature,  the  initial  conditions  for  the  strong-form  problem  (77)-(79)  can  be  stated 
as 


p(x,  0)  =p0(x), 
u(x,  0)  =u0(x), 

0(x,O)  =0o(x), 

for  x  G  12.  The  constitutive  relations  have  been  given  in  Section  2.2.3.  For  the  completeness 
of  this  section,  we  list  them  here: 

r  =p  (Vu  +  VuT)  +  AV  •  ul, 

=  ^ApAp  +  Vp|2^j  I  -  AVp 0  Vp, 

p  =Rb6 — ^ - ap2, 

b  —  p 

q  =  —  kV6, 
n  =ApV  •  uVp. 
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Remark  8.  In  the  remainder  of  this  work,  the  Stokes  hypothesis  is  adopted,  i.e. 


A  =  -3* 

The  total  energy  can  be  represented  as 

X  1 

pE  =  pi  +  -p|u|J  —  +  P®s  +  2dlu|2-  (80) 

The  definitions  of  the  thermodynamic  state  variables  are  recollected  here.  The  Helmholtz 
free  energy  density  T,  the  local  Helmholtz  free  energy  density  \E boo  the  local  internal  energy 
density  iioc,  the  internal  energy  density  1,  the  entropy  density  s,  and  the  local  electrochemical 


potential  zyoc  are  defined  as 

*(/>,0,W>)=«Mp,0)  +  A|Vp|2,  (81) 

2  p 

d' ioc(p ,  0)  =  -  ap  +  R0  log  -  Cv® log  (^~)  +  cv0,  (82) 

L  =^oc+ —\S7  p\2,  (83) 

t'loc  —  ~~  op  +  Cv6 ,  (84) 

S  =  ~  Rlog  +  C°hg  (£)  ’  (®) 
Uloc  =  -  2  ap  +  R9  log  ^ +  ^~~p  -  Cv® log  +  Cv®'  (86) 


3.2  Dimensionless  form  of  the  Navier-Stokes-Korteweg  equations 

In  this  section,  we  perform  dimensional  analysis  of  the  Navier-Stokes-Korteweg  equations 
using  the  MLTQ  system.  The  reference  scale  of  mass,  length,  time,  and  temperature  are 
denoted  as  M0,  L0,  T0,  and  90.  We  may  obtain  the  dimensionless  quantities  denoted  with  a 
superscript  *: 


x  =  Lqx *,  t  =  T0t *,  p  =  ^ p *,  9  =  909*,  u  =  , 

Lo  1o 

_ Mo  x  _  Lq  XHe  _  M0  _* _ M0  ^ 

P  ~  T  7^2  P  ’  A  —  Jl  ,r  rp2  A  1  P  —  T  rp  P  1  T  _  rji2  T  T  > 
MUo  0  MHO  ioivo 

Mo  *  ^ou*  M0L0  Ll 

^  =  ’  b  =  -2b,  k=——k,  E=-^E, 

10L0  Jo  Jo 

Mo  *  -Mo  *  Tg  *  Lg  * 

q  =  ^q  ,  n  =  -n  ,  r  =  -^r  ,  s  =  • 
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With  the  above  dimensionless  variables,  the  dimensionless  balance  equations  can  be  written 


M0  fd(p*u* 


TqLq  V  dt* 
M0  fd(p*E* 


T'nLn  V  3t 


+  v*  •  (p*u*  ®  U*)  +  vy  -  v*  •  r*  -  v*  •  y  -  p*b*J  =  o, 
■  +  V*  •  ((p*£*  +  p*)u*  -  (r*  +  Ou*)  +  v*  ■  q*  +  V*  •  IT 


p*b*u*  -  p*r*  =  0. 


The  constitutive  relations  can  be  rescaled  as 

*  nuL0i*e0Pre*  m0t02  *2 

p  =  Rb-  -  , — — - a — t^c~P  i 

bLl  -  M0p*  Ll 


t*  =  p*  (V*u*  +  V*u*T)  -  -/2*V*  •  u*I, 

o 

y  =  -A *V*p*  ®  V*p*  +  (\*p*A*p*  +  y  |V*p*|2^)  I. 


q*  =  -/c*V*0*, 

i?T20Ol  /  M0p*  ^  ,  CJ^o,  ^0o0^ 

s  =  -^r log  J + -**- log  \<Q )  ’ 

rr  =  a *p*v*  ■  u*v*p*. 

The  dimensionless  viscosity  coefficient  p*  =  LqTqp/Mq  measures  the  ratio  of  the  viscous  force 
to  the  inertial  force;  the  dimensionless  capillarity  coefficient  A*  =  MqTqX/ L‘q  measures  the 
ratio  of  the  surface  tension  to  the  inertia  force.  Hence,  the  two  coefficients  can  be  represented 
in  terms  of  the  Reynolds  number  Re  and  the  Weber  number  We  as 

1  1 

^  ““  Re’  ““  We' 

There  is  one  standard  relation  in  thermodynamics  relating  the  specific  heat  capacity  at 
constant  volume  Cv  and  the  specific  gas  constant  R : 


Cv  = 


7-1’ 
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wherein  7  is  the  heat  capacity  ratio.  Hence,  we  have 


Cv  _  1 

R  7  —  1 

Remark  9.  The  value  of  7  is  related  to  the  degrees-of- freedom  of  the  gas  molecule.  For 
example,  7  for  water  vapor  is  1.33  [62]. 

If  the  reference  scales  are  chosen  as 


Mo  = 

LI  " 

M0 

LqTq 

0o  = 


b, 

ab 2, 


0 


crit 


8  ab 
27 R’ 


and  the  reference  temperature  is  selected  as 


0ref  0criti 


the  dimensionless  Navier-Stokes-Korteweg  equations  can  be  written  as 
SF  +  v*>*u*)  =  o, 


wherein, 


<9(p*u* 

dt* 

d(p*E* 


dt* 


+  V*  •  (p*u*  ®  u*)  +  V*p*  -  V*  •  t*  -  V*  •  <;*  -  p*b*  =  0, 
■  +  V*  •  (( p*E *  +  p*)u*  -  (t*  +  +  V*  •  q*  +  V*  •  II* 


p*b*u*  -  p*r*  =  0, 


P  = 


T  = 


<?  = 


rP* 


27(1  —  p*) 

-J-  (  V*u*  +  V*u*T  -  l V*  •  u*I  )  , 

Re  \  3 

2^  ( (p*A V*  +  1|V>-|2 )  1  -  VV*  0  v>*  ] , 


q*  =  -K*V*0*, 


n.  =  _p.v..u.vv, 


(88) 

(89) 

(90) 


(91) 

(92) 

(93) 

(94) 

(95) 
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(96) 


Re  = 


L0b\/ab 


We=^. 

A 

Likewise,  the  thermodynamic  state  variables  (81)-(86)  can  be  rescaled  as 


^  =  ^0C(P^)  +  4e^|VV|2, 


yUp,e)  =  -p*  +  -6*\og 


1  -P*J  27(7-1) 


e*  log (9*)  + 


27(7  -  1) 


<*  =  llr  + 


^W|Vyf’ 


*  _  *  |  /l* 

hoc  -  -P  +  27(7  —  1)  ’ 

*  _  o  *  ,  80*  ,  8  m  i  (  P * 

Uloc  P  +  27(1  -  p*)  +  276  °g  \1  -  p* 


27(^I)ri°g^  + 27(^1)  ^ 


r,  (99) 
(100) 
(101) 


s*  = - log  — -  +  — - log  (0*) . 

27  h\l-p*J  27(7-1)  ; 

Remark  10.  The  dimensionless  thermal  conductivity  k*  can  be  rewritten  as 


(102) 


(103) 


*  90Tq  8  n  8  k  8  7  1  k  8  7  1 

MqLq  K  27  RL0bVab  27  RRefi  27  7  -  1  Re  Cpp  27  7  -  1  Re  Pr  ’ 

wherein  Cp  :=  7 Cv  is  the  specific  heat  capacity  at  constant  pressure  and  Pr  :=  Cppj k  is  the 
Prandtl  number  [56]. 

Remark  11.  The  capillarity  number  Ca,  which  measures  the  relative  effect  of  the  viscous 
force  against  the  surface  tension,  is  defined  as 

We 

Ca  =  l- 

The  Bond  number  Bo  measures  the  ratio  of  the  body  force  to  the  surface  tension  and  it  is 
defined  as 

Bo  =  lb* I  We . 


7  1 


Henceforth,  we  will  restrict  our  discussions  to  the  dimensionless  form,  and  the  superscript 
*  will  be  omitted  for  notational  simplicity. 
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3.3  Functional  entropy  variables 

The  mathematical  entropy  function  H  is  defined  to  be 


(104) 


With  this  definition  and  the  relation  (72),  the  second  law  of  thermodynamics  can  be  written 
in  terms  of  H  as 


f +  V.(i4»)-V.(|)  +  ^  =  -ir:Vu-iKW<0. 


In  three  dimensions,  the  conservation  variables  can  be  written  as 


UT=  [U1:U2,U3,U4:U5\  :=  [p,pu1,pu2,pu3,pE]. 


The  classical  entropy  variables  for  the  compressible  Navier-Stokes  equations  are  defined  as  the 
partial  derivatives  of  the  mathematical  entropy  function  H  with  respect  to  the  conservation 
variables  U.  This  definition  of  the  entropy  variables  was  understood  as  an  algebraic  change- 
of-variablcs,  since  the  mathematical  entropy  function  for  the  compressible  Navier-Stokes 
equation  is  a  function  of  the  conservation  variables.  In  contrast,  due  to  the  constitutive 
relation  (100),  the  temperature  9  for  the  van  der  Waals  model  can  be  expressed  in  terms  of 
the  conservation  variables  as 


9  = 


27(7  -  1)  f*h 

Fi 


uj  +  Uj  +  Uj 

2  U\ 


2  Webb 


\VU1\2  +  U1 


The  above  relation  includes  a  non-local  gradient-squared  term.  This  fact  suggests  that  when 
taking  derivatives  of  the  temperature  with  respect  to  conservation  variables,  the  deriva¬ 
tion  should  be  taken  in  the  functional  setting.  Therefore,  for  the  Navier- Stokes- Korteweg 
equations,  we  define  the  entropy  variables  V  as  the  functional  derivatives: 


V 


5H 

hU 


Wi'i  V2;  V3;  14;  V5]T 


'  5H  5H  5H  5H  5H 
Wi  W2  W3  Wi  6U~5 


T 


Given  the  test  functions  <5v  =  [<5ui;  Sv2  \  Sv3,  Sv 4;  5v5]T,  the  entropy  variables  V  can  be  written 
explicitly  as 


4,wd(-2p+^oe(^ 


1  ~pj  27(7-1) 


9  log  (9)  + 


27(7  -  1) 


9 
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8  9  U2\  £•  1  1,-,  nt 

+  27(l-p)  2  + 

(105) 

V2[Sv  2]  =  y<k>2, 

(106) 

V3[Sv3]  =  ^-8v3, 

(107) 

Va[5va\  =  y^’4, 

(108) 

V5[Sv5]  =  ~^V5- 

(109) 

Remark  12.  The  local  Helmholtz  free  energy  p40oc  can  be  regarded  as  a  function  of  p  and 
9.  Taking  derivatives  of  p^>ioc  gives 


H 


L'loc 


d(p^ioc) 

89 

d(p^ioc) 

dp 


Remark  13.  In  Section  2-4-2,  we  introduced  the  local  electrochemical  potential  vioc,  defined 
as 


d(p^ioc)  dtyioc 

Vioc  =  - a -  =  ^ loc  +  P~ 


dp 


dp 


We  define  the  global  electrochemical  potential  v  by  generalizing  the  partial  derivative  in  the 
above  formula  to  the  functional  derivative: 


5^ 

v[Sv i]  :=  ^[Svi]  +  p—[Sv i] 

op 


=  (  -2 P+  log 


P 


1  -P)  27(7-1) 


6  log  (9)  + 


27(7  -  1) 


9 


27(1 -p) 


dvi  +  — Vp  •  V5v i, 
We 


=  vloc5v i  +  —  Vp  •  V6v i. 

We 

hiterestingly,  with  this  definition,  the  entropy  variable  V\  can  be  written  as 


Vi[Svi\  =  g  (  v 


u 


M- 
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Consequently,  the  entropy  variables  V  can  be  compactly  represented  as 


v  — 


V  = 


1 

6 


u  l 
u2 


«3 
-  1 


(110) 


The  expression  (110)  formally  coincides  with  the  definition  of  the  entropy  variables  for  the 
perfect  gas  model.  However,  the  entropy  variables  here  should  be  understood  as  linear  oper¬ 
ators  in  the  dual  spaces  of  the  conservation  variables.  The  expression  (110)  also  hints  that 
the  formulation  of  the  entropy  variables  is  invariant  under  different  choices  of  the  Helmholtz 
free  energy  functional. 


Theorem  3.  The  action  of  entropy  variables  V  on  the  Navier- Stokes- Korteweg  equations 
recovers  the  Clausius- Duhem  inequality. 


Proof.  Testing  the  entropy  variables  V  with  the  time  derivative  terms  leads  to 


~d\J~ 

5H 

'0U' 

.  dt . 

'  6V 

.  dt  . 

d  H 
~dt 


Choosing  the  test  function  as  the  Euler  flux  results  in 


(111) 


V 


V  •  (pu) 

V  •  (pu  (g)  u)  +  Vp 

V  •  (pEu  +  pu) 


=  V  ■  {Hu) 


Wed 


Vu  :  Vp  <g>  Vp 


Vp|2V  •  u  +  pVp  •  V(V  •  u) 


(112) 


We  emphasize  that  the  notation  V  [  ]  denotes  the  action  of  the  differential  operator  V  on 
the  term  in  square  brackets.  Taking  the  test  functions  as  the  terms  related  to  the  capillarity 
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leads  to 


—V  •  fcu)  +  v  •  n 


i  /  i  _ 

XYvll  (  Vu:  Vp<g>Vp  +  -|Vp|2V-u  +  pVp-  V(V-u) 


(113) 


Combing  (112)-(113)  yields 


V  •  (pu) 

V  V  ■  (pu  ®  u)  +  Vp  —  V  •  <?  =  V  ■  (Hu) 

V  ■  ( pEu  +  pu)  —  V  •  (<?u)  +  V  •  n 


(114) 


Testing  the  entropy  variables  against  the  viscous  flux  gives 


V  -  V-r 
—  V  •  (ru) 


=  r  ■ Vu- 


(115) 


The  action  of  entropy  variables  on  the  heat  flux,  the  heat  source,  and  the  body  force  yields 


—  pb 

V  ■  q— pu  •  b  —  pr 


=  -v'(?)  +  f+^'w  I2- 


(lie) 


Combing  the  relations  (111),  (114),  (115),  and  (116)  leads  to 


f +  V.(»u)_V.(5)+^  =  ,^|Vfl|»-lr:Vu, 


or  equivalently, 


+  v  .  („su)  +  V  ■  (3)  -  £  =  IK|V0|2  +  \r  :  vu. 
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This  is  exactly  the  dissipation  relation  for  the  Navier-Stokes-Korteweg  equations.  □ 


3.4  An  alternative  statement  of  the  strong  form  of  the  problem 

Theorem  3  suggests  that  a  weak  formulation  for  the  Navier-Stokes-Korteweg  equations  will 
satisfy  the  Clausius-Duhem  inequality  weakly  as  long  as  the  entropy  variables  V  are  in  the 
test  function  spaces.  For  the  compressible  Navier-Stokes  equations,  one  may  rewrite  the 
equation  in  terms  of  the  entropy  variables  V,  since  the  mapping  between  U  and  V  is  purely 
algebraic.  By  using  the  Bubnov-Galerkin  method,  the  entropy  variables  are  enforced  in 
the  test  function  spaces,  and  consequently  one  can  prove  the  entropy  stability  for  the  finite 
element  formulation  [34],  This  approach  has  been  adopted  for  constructing  entropy  stable 
finite  element  formulations  for  a  variety  of  problems  [4,  8,  10,  34,  35].  However,  for  the 
van  der  Waals  fluid,  there  is  an  additional  difficulty  coming  from  the  differential  relation 
in  the  definition  of  V\  in  (105).  The  classical  approach  becomes  nonviable,  since  there  is  a 
second-order  differential  operator  in  the  definition  of  Vf,  and  inverting  a  differential  operator 
is  not  a  straightforward  task.  Inspired  from  the  form  of  V) ,  we  introduce  a  new  independent 
variable  and  couple  it  with  the  balance  equations  by  replacing  the  pressure.  Hence,  we  may 
derive  a  new  system  of  equations,  which  is  a  consistent  statement  of  the  original  strong- 
form  problem.  In  doing  so,  the  definition  of  the  entropy  variable  V\  is  weakly  enforced  for 
the  mass  balance  equation,  and  we  can  prove  entropy  stability  for  the  corresponding  weak 
formulation.  To  derive  the  alternative  statement  of  the  Navier-Stokes-Korteweg  equations, 
we  introduce  the  auxiliary  variable  V  as 


Recall  that  the  local  electrochemical  potential  is  related  to  the  thermodynamic  pressure  by 


p 

L'loc  1“  ® loC' 

p 

Hence,  the  auxiliary  variable  V  can  be  rewritten  as 


V 


1 

0 


+  ^loc  — 


Rearranging  terms  in  the  above  relation  yields 


P  =  pVd  -  p^ioc  + 


p|u|2 

2 


+ 


1 

We 


p0X7  ■ 


(118) 
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The  above  relation  is  an  equivalent  expression  of  the  van  der  Waals  equation  of  state  (67)  in 
terms  of  the  newly  introduced  auxiliary  variable  V.  Taking  gradient  of  both  sides  of  (118), 
we  have 


V  (pdhoc) 

d{p^ioc)  _  8  (pTfOC) 

dp  p  89 


=  v  (pVB  +  -  HV0 


(119) 


Using  (118),  the  term  pE  +  p  can  be  reorganized  as 

pE  +  p=p^ioc  —  9H  +  ^-|Vp|2  +  ^p|u|2 
+  pVO  +  -p|u|2  —  P^ioc  +  ^-p^V  ■ 

-  SH  +  s4e  I ’ +  ^lul2  +  '  (t)  ‘  (120) 


Making  use  of  (119)  and  (120),  the  pressure  force  Vp  and  the  power  expenditure  of  pressure 
V  •  (pu)  can  be  consistently  represented  in  terms  of  V.  Then  the  original  strong-form  problem 
(88)- (90)  can  be  rewritten  as 


(121) 

(122) 

(123) 

(124) 
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The  equation  (124)  defines  the  auxiliary  variable  V.  Based  on  our  derivation,  the  new 
balance  equations  (121)-(123),  together  with  the  auxiliary  variable  (124),  is  equivalent  to 
the  original  Navier-Stokes-Korteweg  equations. 

3.5  Weak  formulation 

In  this  section,  we  construct  a  weak  formulation  based  on  the  alternative  statement  (121)- 
(124).  In  the  weak  formulation,  we  solve  for  six  unknowns  in  three  dimensions.  The  set  of 
variables  is  denoted  as 

P 

Ui 

e 

U2 

e 

U3 

e 
i 
e 

v 

Let  Vi  be  the  trial  solution  space  for  Y\  —  p  and  Y§  —  V;  V2  be  the  trial  solution  space  for 
Yi+ 1  =  Ui/9,  i  =  1,  2, 3;  V3  be  the  trial  solution  space  for  Y'5  =  —1/6.  The  test  function  spaces 
are  taken  to  be  identical  to  the  corresponding  trial  solution  spaces.  The  weak  formulation  can 
be  stated  as  follows.  FindYi(t)  =  pit)  G  L2(0,  T;  Vi)fli/1(0,  T;  L2(fl)),  Yi+i(t)  =  Ui(t)/9it)  G 
L2(0,  T;  V2)nH1iO,  T ;  L2( Q))  for  i  =  1,  2, 3,  Ys(t)  =  -1  /9(t)  G  L2(0,  T;  V3)nH\0,  T;  L2( ft)), 
and  Y6(t)  =  V  G  L2{ 0,  T;  Vi),  such  that 

wi,  ~  {Vwup u)n  =  0,  Vwi  G  Vi,  (126) 

(w-^ir)n- (Vw- pu  ®  u)n "  (v ' w- pVe + y |u|2  +  ■  (t)  )  n 

-  (w,  (v9  +  Vp^)  -  (w ,HV9)n  +  (Vw,r)n  +  (Vw,?)fl 

=  (w,pb)n,  Vw  =  (1V2',  w/s/)T  G  (V2)3  ,  (127) 
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(W5’^)n"(Vl"5’H-m+^|Vp|2+',|u|2+^v'  (t)) 

+  (Vw5,  rujjj  +  (Vic5,  <?u)n  -  (Vw5,  q)n  -  (Vw5,  n)n 
=  (w5,  pb  •  u)n  +  (w5,  pr)n  ,  Vu>5  G  V3, 

(”’»•  =  (“«’  l  -  ^-)  )  „  +  (Vw‘’  vLV/>)n’  V!"s  6  Vl' 


n 


(128) 

(129) 


with  p(0)  =  po,  u(O)/0(O)  =  uo/^o,  and  —1/0(0)  =  —  l/0o  in  fh 

Comparing  (129)  with  (105),  one  may  find  that  the  auxiliary  variable  V  is  identical  to  the 
entropy  variable  V\  in  the  weak  formulation.  Therefore,  in  the  set  of  variables  (125),  we  are 
actually  solving  for  the  entropy  variables  V  together  with  density  p,  which  is  the  conjugate 
variable  to  Vj  =  V.  By  choosing  the  test  function  and  trial  solution  spaces  identical  for 
the  equations  (126)  and  (129),  the  entropy  variable  V\  is  weakly  enforced  to  be  in  the  test 
function  space  for  the  mass  conservation  equation.  This  is  a  key  ingredient  in  the  proof  of 
the  following  theorem. 


Theorem  4.  Sufficiently  smooth  weak  solutions  of  the  problem  (126)-(129)  verify  the  second 
law  of  thermodynamics,  i.e., 


_f  +  v.„uKV.(3)  +  f|^ 


/  -t  :  VudVx 
In  0 


K,|V0/ 

~lP~ 


dVx.  (130) 


Proof.  Choosing  w\  —  V  in  (126)  and  w6  =  dp/dt  in  (129)  yields 


dp  y\  _  (dp  1  /  _  Jup 

Of'  )v  {ot'e{Uloc  2 


+  I  V  |  .  f-^P 


dt  J  ’  W eO 


Combing  the  above  two  relations  leads  to 
dp 


r  sh 

In  dp 


dt 


d\f-  [uloc 

—  (VV,  pu)n . 


u 


V|3? 


dt  J  '  Wed 


Vp 


(131) 


Taking  w  =  u/6  in  (127)  results  in 


d(pu) 

dt 


d\ 4 


u  <9(pu)  \ 
0’  dt  )n 


(V(|),pu®u)n 
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+  2  +  G^V'(t 


u  /„„  iur  l _ (S/p 


«■  M+y+#'  t'IVp 


n 


?-ffV9)fi-(v(?)'")n-(vG)'Sn 

^b)v 


Choosing  w5  =  —1/6  in  (128)  yields 


r  5h  r  d{pE) 

In  Kpe)  .  dt 


=-(v(l],i  pvs-oh 


e  ’  at  y 


+  24e|V',|2  +  p|u|2  +  ye',9V'(irl  1“ 

vGWfi+(v6Kn 

5,v‘q)0"(vG)'n)0"G'',b'u/o 


0’ 


pr 


Grouping  all  terms  in  (131)-(133)  involving  V,  one  has 

(W,pu)„ 


vG)'Hn+G’l/Hn-(vui'l/^ 


=  /  V-(pVu)dVx  =  /  pVu-ndAx  =  0. 
Jn  Jan 


Summing  all  terms  in  (131)-(133)  involving  H  yields 

(?'flv0„+(vG)'9j/u)n= g-hw)«-  G’v(wu) 

=  -  [  V  •  (Hu)  dVx. 

Jn 

Next,  collecting  all  the  terms  in  (131)-(133)  explicitly  involving  We,  we  have 


-O'®”' 

-  (  v  f-']  —  P6V  ■  (— \  Y  -Yvf-V  1  IV7-'2’ 


n 


9  I  We 


6 


6  7  2  We 


I Vp|  u 


(132) 


(133) 


(134) 


(135) 
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=  /  ^-V  ■  (pu)  v  ■  (  V  ]  -  —  v  U  ]  ■  u|Vp|W; 


We 


e 


2  We  \d 
1 


Wed 


-1 
We9  \  2 


V  (V  •  (pu))  •  Vp  +  ^-^V  ■  (u|Vp|2)  dVx 
|Vp|2V  ■  u  +  Vp  <g)  Vp  :  Vu  +  pVp  •  V  (V  •  u)  )  dVx. 


(136) 


Combing  all  the  terms  in  (131)-(133)  including  the  Korteweg  stress  and  the  interstitial 
working,  we  have 


___  i  u  , 


V  (  £  ) 


vU)'n 


d - 


Vu  :  I  +  -,V-n 


d 


-|  Vp|2V  ■  u  +  Vp  ®  Vp  :  Vu  +  pVp  •  V  (V  •  u)  )  dVx. 


Jn^ed  \2 

Making  use  of  (134)-(137),  the  summation  of  (131)-(133)  gives 

dVx 


r  sh 

dp 

5H 

1 

'(9(pu)' 

5H 

\S(pE)  1 

L  sp 

_dt_ 

5(pu) 

dt 

8{pE) 

dt 

(  V  •  (Hu)  dVx  —  (V  (^),- 
J  n  \  \t>  / 


-,pb-u  +  -,V-q 


7T  »T  _+  V  L  -TU  +  -,pb-u 


d 1 


n 


d ’ 


0’ 


0’ 


pr 


n 


=  -  I  V  •  (i/u)  dVx  -  f  Vu dVy 

Jn  Jn  & 

=  -  I  V  •  (i/u)  d\ 4  ~  I  VudVx  +  /  V 
Jo  Jo  d  Jq 

The  above  equation  implies 


+ '  rv'q 


q 


0’ 


pr 


o 


?l  + 


q-  _  P^ry 

d2  d  x' 


(137) 


which  completes  the  proof  of  this  theorem.  □ 

Remark  14.  In  our  discussion,  we  assumed  periodic  boundary  conditions.  For  periodic 
boundary  conditions,  the  divergence  terms  in  (130)  are  canceled  out  and  the  statement  can 
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be  simplified  as 


1 

9 


t  :  \7udVx  — 


n\V9\2 

92 


d\ 4- 


Even  though  we  proved  the  case  with  periodic  boundary  conditions,  the  proof  of  Theorem 
4  can  proceed  with  other  boundary  conditions,  such  as  no-slip  boundary  conditions  for  the 
velocity  field  and  the  heat  flux  boundary  condition  for  the  temperature  field. 


3.6  Semi-discrete  formulation 


We  perform  spatial  discretization  of  (121)-(124)  by  invoking  the  Galerkin  method  [32],  Let 
Vf  C  Vi,  V2  C  V2,  and  V3  C  V3  be  finite-dimensional  function  spaces,  in  which  the  super¬ 
script  h  denotes  a  mesh  parameter.  Then  the  spatial  discretization  of  (121)-(124)  can  be 
stated  as  follows. 

Find  Yfiit)  =  ph(t)  G  L2(0,T-Vf)nH\0,T-L2(n)),  Yf+1(t)  =  u^{t)/9h{t)  G  L2(0,T;V£)  n 
H\0,T-,L2(n))  for  i  =  1,2,3,  Y5h(t )  =  -l/9h(t)  G  L2(0,  T;  )  D  H\0,  T;  L2(fi)),  and 
Yfi  =  Vh  G  L2(0,  T;  Vf ),  such  that 


w 


w 


f,  -  (Vw},puh)a  =  0.  6  V*. 


d(phu.h) 

dt 


n 


V  •  wh,  phVh9h  +  l-ph\nh\2  +  ^rPh9hV  ■  ( ^ 
2  We  \  9n 


,h  12 


w'*.  ( ^ 


(w  h,HhV9h)i 

+  (Vwfc,Th)n+  (Vwk,sk)n  =  (w\phb)n,  Vw''  =  (w%;w!l;w!t)T  G  (V2fe)3  , 
h  _  (Vu;5’  - ehHh  +  +  ph|uT 


w 


5> 


1  \7  f 


uM  +  (VW5,  +  (VWg, 


-  (Vwi,  q'%  -  (Vwi,  n'%  =  (w»  p%  ■  u'%  +  (w‘,  A)„ ,  M  6  Vj, 


W,  On  =  («&  f  ^ 


+  |VW«’We^ 


w 


M  6  vf, 


(138) 


(139) 


(140) 

(141) 


43 


wherein 


h - yh 

P  ■—  1 1  i 


uh  :  = 


6h  ■= - -r, 

Vh  ’ 


(  y-t 

Vh 
I3  . 

V  h 
1 4 

V  YP' 

Yf  ’ 

Yf 

r/!:=R^(Vu/l  +  (Vu^T“^v'u"1)’ 


phAph  +  ^\Wph\2  I  -Vpft®Vp‘ 


h  1 2 


We 

qh  :=  -KV0h, 

nft  :=  -*-pftV  •  uhX7ph, 
We 

ph 


^  :=  ^  (  log 


7-1 


<,:=-2ph  + 


1  —  ph 
86h  8  nh  , 
27(1 -ph)  +  27^  °S 

9h  \og(dh)  + 


log(^ 


P 


27(7  -  1) 


27(7  -  1) 


P' 

0fc, 


with  ph( 0)  =  Pq,  uh(0)/6|/l(0)  =  Uq/Pq,  and  —  l/6*h(0)  =  —  1/0q  in  fl. 

In  the  above  formulation,  pp,  Uq/Pq,  and  —1/0q  are  L2-projections  of  p0(x),  uo(x)/0o(x), 


m 


and  — 1/6*0  (x)  onto  Vy,  'Ey,  and  ’Ey  respectively.  Employing  the  same  techniques  used  i 
the  proof  of  Theorem  4,  we  can  obtain  the  following  theorem,  which  implies  that  the  spatial 
discretization  (138)-(141)  is  entropy  dissipative. 

Theorem  5.  The  solutions  of  the  semi-discrete  formulation  (138) -(141)  satisfy  the  second 
law  of  thermodynamics  in  the  following  sense. 


Remark  15.  In  our  implementation,  the  same  discrete  space  Vh ,  up  to  the  prescription 
of  the  boundary  conditions,  is  used  to  approximate  V\ ,  V2,  and  V3.  Specifically,  the  Non- 
Uniform  Rational  B-Spline  (NURBS)  basis  functions  are  used  to  define  Vh  as  well  as  the 
geometry  of  the  computational  domain.  Consequently,  this  approach  may  be  considered  as 


44 


isogeometric  analysis  method  [33]. 


3.7  The  fully  discrete  formulation 

In  the  preceding  section,  we  have  constructed  an  entropy-dissipative  semi-discrete  formula¬ 
tion.  It  remains  to  design  discretizations  of  the  time  derivatives  such  that  the  dissipation 
property  can  be  inherited  in  the  time  direction.  In  our  previous  work  [45],  we  have  suc¬ 
cessfully  developed  a  suite  of  temporal  schemes  for  the  isothermal  Navier- Stokes- Korteweg 
equations.  For  the  thermal  case,  the  difficulty  comes  from  the  term  d(pE)/dt.  If  one  uses 
the  traditional  jump  operator  to  approximate  the  time  derivative,  it  is  difficult  to  estimate 
the  dissipation  of  the  resulting  scheme.  In  this  work,  the  total  energy  pE  is  split  into  four 
parts: 

p£  =  ^te-«/+ip|u|2+^_|V>|2,  (142) 

and  the  time  approximation  for  each  of  the  four  parts  will  be  carefully  designed  to  ensure 
consistency  and  temporal  dissipation.  It  is  noteworthy  that,  in  the  design  of  the  discrete 
scheme,  the  special  quadrature  rules  developed  in  [24,  45]  will  be  used  repeatedly  as  a  key 
technique.  In  the  following,  we  will  first  state  the  fully  discrete  scheme  in  Section  3.7.1. 
Following  that,  five  preliminary  lemmas  will  be  given  in  Section  3.7.2.  The  main  results 
about  the  entropy-dissipation  property  and  time  accuracy  are  proven  in  Section  3.7.3. 


3.7.1  The  fully  discrete  scheme 

To  discretize  the  semi-discrete  formulation,  the  time  interval  X  =  (0,  T)  is  divided  into  Nts 
subintervals  Tn  =  (tn,  tn+ 1),  n  —  0,  •  •  •  ,  Nts  —  1,  of  size  A tn  =  tn+ 1  —  tn.  We  use  the  notation 


-yri 
x  n  * 


h.  Ul  n. 

Pm  nh  ’ 

wn 


«2  ,n. 

eh : 

n 


U 


Ah-  _L  y 

eh  ’  eh  ’  * 

n  n 


(143) 


to  denote  the  fully  discrete  solutions  at  the  time  level  n.  The  fully  discrete  primitive  variables 
at  the  same  time  level  can  be  represented  in  terms  of  Y[]  as 


Phn  =  Ph(Yhn )  =  Y\hn, 

<n  =  <(Yhn)  =  *Y*1>n/Y*n,  i  =  1,2,  3, 
On  =  0h(Yhn )  =  -1  /Y*n. 
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We  define  the  jump  of  density,  linear  momentum,  and  total  energy  over  each  time  step  as 


[Pnl  :=Pn+ 1  -  Pn , 
lPnUn  1  -=Pn+ lUn+l  ~  PnUn, 

[phnE(phn,uhn,9hn)]  :=(P*ioc)(phn+i,Ohn+ 1)  -  (p^loc)(phn+h,Ohn) 

+  (p*ioc)(phn+1,ehn+h)  -  (p*ioM,ehnH) 
-0hnH  ( H(phn+l,ehn+1)-H(phn,9hn )) 

_  Ci_c:  (ff(p;+,,»;+I)+jf(P;+,,o;)) 


'  \P n_y_  1  5  ^71+1/ 


12 


+  ^(p: 

1 


n+l  Wn+1 


1 2 


2  We 


(|Vp: 


h  | 
n+l  I 


p>J!2) 


VC) . 


(144) 

(145) 


(146) 


Remark  16.  According  to  the  energy  split  (142),  the  definition  (146)  can  be  rewritten  as  a 
summation  of  four  jumps: 


wmpI  <  o]  «;;)]  -  oi  +  i^i<ri  +  i^r  r vci,  (147) 


wherein, 


\pyupier,] :=(p«,M)(p» 


(^loaKA-l.^+l)  -  (P«foc)(CCi 


KfCDdi  :=0,';+i  (idco  - 

'  $hn+l~e"  (H(p^,,ehn+1)  +  H(PhnH]e':)) 


2 


r+ 
1  2 


u 


/i  |2i 


12 

1  + 

n+l  I  un+l 


89 2 


n+ 


.Cl), 


I2  —  pt|uAl2 


), 


IT— i— IVdl  : 
l2We  rA  * 


2  We 


(IVpCl2  -  VC) 


(148) 


(149) 

(150) 

(151) 


The  definitions  (148)  and  (149)  are  inspired  by  the  fact  that  the  summation  of  first-order 
partial  derivatives  approximates  the  total  derivative.  In  (149),  there  is  an  additional  third- 
order  perturbation  term,  whose  role  will  be  described  in  Lemma  6.  The  jump  operators  for 
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the  kinetic  energy  and  the  surface  energy  follow  the  classical  definition.  In  the  following  text, 
we  define  :=  9^+1  —  9 to  simplify  notations. 

With  the  jump  operators  defined  above,  the  fully  discrete  scheme  is  stated  as  follows. 
In  each  time  step,  given  and  the  time  step  A tn,  find  Y(j+1  such  that  for  all  G  Vh, 
wh  =  (w%]  w 3;  wf)T  G  (Vft)3,  wQ  G  Vh,  and  G  Vh, 


.h  iPn  1 


BmK;Y^+1)  :=  (<  At 
B['(w*;Y;+1)  :=  (V, 


n  /  n 
Jh  I  PnUnl 


-  (vwf.pJLi , 


n+2  n+2  /n 
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_2_  I  Qh  V7  I  lt+2 
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Vp 
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+  (Vw  ft,TA 
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n+Uo 


0 

-  (wft,Pn+ib)^  =  0,  (153) 
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[PnF^oci  ft  fth 
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12  dp 2 
wherein 


(152) 


(154) 


=  0,  (155) 


Yn+| 

dn+I  - 

uh  1 

n+j 


2  (Yn  +  Yn+l)  > 

p\  y;+i), 


:=  u^fY71 


n+1 


n+5  ’ 


=  ^(Yn+l), 


(156) 

(157) 

(158) 

(159) 
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(160) 


(161) 

(162) 

(163) 

(164) 


3.7.2  Preliminary  lemmas 

We  state  and  prove  five  lemmas  in  this  section,  which  will  be  applied  to  prove  the  final 
results  in  Section  3.7.3. 


Lemma  3.  The  mathematical  entropy  function  H(p,  6)  given  by  (104)  satisfies 


d3H 

~w 


<  0. 


Proof.  Straightforward  calculations  lead  to 

d3H  16  p 

~d93  ~  _ 27(7  —  1)  d3' 


(165) 


The  dimensionless  temperature  9  is  always  positive  and  the  heat  capacity  ratio  7  is  always 
greater  than  1.  Hence, 


d3H 

~W 


<  0. 


□ 


Lemma  4.  The  local  electrochemical  potential  z7oc(P)$)  given  by  (102)  satisfies 


Proof.  Direct  calculations  yield 


d3Poc 

dp3 


>  0. 


d3vioc  16 9  6 p2  —  Ap  +  1 

dp3  ~  ~W  p3(l-p)4 


(166) 
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It  is  known  that  6  >  0  and  6 p2  —  4p  +  1  >  1/3.  Therefore,  one  has 


d\oc 

dp3 


>  0. 


□ 


Lemma  5.  (Perturbed  trapezoidal  rules)  For  a  function  f  G  C3([m,n]),  where  m,n  €  M, 
there  exist  £i,£2  G  ( m,n )  such  that  the  following  quadrature  formulas  hold  true. 


ft  \j  n~m  (ft  \  ,  ft  w  in-mf  f"(  (n  -  m)4 
f(x)dx  =  — —  (f(m)  +  f{n )) - — - f  (m) - — - /  (£1), 


12 


24 


77  —  777 

f(x)dx  =  — —  (f(m)  +  f(n)) 


(n-m)3  „  (n-m)A  „, 

f  (n)  + - 55 - f  (&)■ 


12 


24 


(167) 

(168) 


The  proof  for  this  lemma  can  be  found  in  the  appendix  of  [24].  There  are  two  other  suites 
of  quadrature  formulas  -  the  rectangular  quadrature  rules  and  the  perturbed  mid-point  rules. 
Interested  readers  are  referred  to  [45,  43]  for  details  about  these  formulas  and  applications.  A 
common  feature  of  these  formulas  is  that  each  pair  contains  opposite  signs  in  the  asymptotic 
residual  terms.  This  allows  one  to  perform  a  split  of  the  target  function  and  construct  a 
discrete  scheme  with  a  controllable  residual.  This  technique  will  be  demonstrated  in  the 
following  lemma. 


Lemma  6.  Given  Ipfj,  [Pnun],  and  [p^E(p^  u£,  0^)]  defined  in  (144)-(146),  the  following 
relation  holds  for  £i,£2  £  (0, 1). 
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(169) 
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Proof.  Direct  calculations  can  show  that 
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Making  use  of  the  above  two  relations,  one  can  get 
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Applying  the  perturbed  trapezoidal  rule  (167)  to 
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one  can  get 
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for  £1  e  (0, 1).  Consequently,  the  relation  (170)  can  be  rewritten  as 
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Next,  applying  the  perturbed  trapezoidal  rule  (168)  to 

rr  d(p^ioc) 


leads  to 
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for  ^2  ^  (0, 1).  Using  the  above  relation,  relation  (171)  can  be  further  rewritten  as 
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9^114  03  TT 
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This  completes  the  proof  of  Lemma  6. 


□ 


Remark  17.  Based  on  Lemmas  3  and  4,  one  can  show  that  the  last  two  terms  in  (169) 
satisfy 


i  Atn’ 
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3^114 
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The  two  terms  represent  the  numerical  dissipation  introduced  by  the  approximation  of  time 
derivatives. 


Remark  18.  In  the  proof  of  Lemma  6,  it  is  clear  that  the  novel  jump  operator  (146)  is 
designed  based  on  the  perturbed  trapezoidal  formulas  (167)  and  (168).  It  should  be  pointed  out 
that  one  may  construct  different  discrete  jump  operators  for  d(pE )  / dt  by  using  the  perturbed 
mid-point  rules  or  the  rectangular  quadrature  rules  proposed  in  [45].  The  resulting  schemes 
can  be  shown  to  guarantee  the  dissipation  property,  but  the  amount  of  numerical  dissipation 
will  be  slightly  different. 

Lemma  7.  Replacing  p vff,  and  91]  in  the  definition  (146)  with  corresponding  time  con¬ 
tinuous  functions  ph{tn),  u h(tn),  and  9h(tn )  and  assuming  sufficient  smoothness  in  the  time 
direction,  one  has 


L ph{tn)E(ph(tn ),  u  h(tn),  eh(tn ))] 

=  (ph(tn+1)E(ph(tn+1),uh(tn+1),Oh(tn+1))  -  ph(tn)E(ph(tn),uh(tn),Oh(tn))^)  +  0{At3n). 

Proof.  Recalling  the  relations  (147)-(151),  we  only  need  to  analyze  the  non-classical  jump 
operators  (148)  and  (149).  We  consider  the  jump  operator  (148)  first.  Taylor  expansions 
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lead  to 
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(p>P1«)(p',(t„+i),9'‘(«„+1))  =  (P'Pi«)(p'‘(i„+i),«'‘(i„+i)) 

+  ^^(P'“(i„+i),«'“(i„+i))  («'‘(«„+1)  -0*(i„+.)) 

+  5^^(p'‘«~+i).',',«n+i))  («"(«,.«)  -«"(i„+i))2 
+  G(At3),  (174) 

(p^oc)(p'l(tn+i),^(tn))  =  (p^oc)(ph(tn+i),^(tn+i)) 

+  ^pV(i„+.  ),«'■«„+.))  («'*(«„)  -  «'■(«„+.)) 

+  -  ^(‘n+l))2 

+  O(Atl),  (175) 

(f,<Uloc)(ph(tn+l),t)h(tn+ .))  =  (p'tI«)(p'‘(i„+.),»',(i„+.)) 
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+  d-^^(ph(tn+l),eh(tn+ ,»  (p‘((„+1)  -  Ph(tn+l)) 

+  0(At,s,),  (176) 

(f>4'te)(/(«„),9',(*„+i))  =  (p't,ocKpb(tn+i),eh(tnH)) 

+  ^^(p'“(i„+i),  «'“(*„+!))  (/(in)  -  /(<„+.)) 

+  ^0^(p'‘(f„+|)-«'‘(‘„+|))  (/W  -  AW l))2 

+  O(AjJ).  (177) 


Combining  the  above  Taylor  expansions  leads  to 

((P®k«)(/(in+l),#*(t«+.))  -  (P*loc)(p'*(i„),«'‘(«„))) 

-  {(p<i,oc)(ph(tnH),eh(tn+1))  -  (pffloc)(ph(tnH),eh(tn)) 

+  (p®i(J(AWi),^«„+i))  -  (p*;oo)(p'‘(f„),  eh(t„+i)))  =  O(Ati). 

Next,  we  analyze  the  term  (149). 


ef‘(tn+1)H(ph(tn+l),el'(tn+1))-eh(tn)H(Pb(tn),eh(tn))') 

-  b*(«„+.) (h (/(«„„), »b(t„+1))  -  H(Ph(tn),eh(tn))) 

h  °h{u+')-  0h{tn)  (n(ph(tnH),eh(tn+1))  +  H(/‘(tn+i),eh(Q)) ) 


eh(tn+1)  -eh(tn)  /  h  h 


(tn+1),e\tn+1))  +  H(Ph(tn),e\tn)) 


-  H(Ph(tn^),eh(tn+1))  -  H(Ph(tn+h),eh(tn)))  +  o(Atl) 


(V(tn+I)  -e\tn))  dHIJtu  w,,. 
- 2 - \-dp(p  (t"+Vie 


Pk(tn+ 1)  -  Ph(tn+  i 


^(/(W‘).^(4))  (A*»)  ~  A^)J  I  +0(A£) 


=(T(A4). 


(178) 
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According  to  (172)  and  (178),  it  can  be  concluded  that 


[ph(tn)E(ph(tn),  u h(tn),  6h(tn))j 

-  ^ oh(tn+l)E(ph(tn+1),uh(tn+1),9h(tn+1 ))  - 

o,  (9h(tn+1)-9h(tn)fd2H/h/ 
=0(At3n)  +  1  Kn+1\n  V  ”  ^{Ph{tn+ 1), 


12 


c>02 


=d(A^ 


ph(tn)E(ph(tn),uh(tn),9h(tn ))) 

0h(Wi)) 


This  completes  the  proof  of  the  lemma.  □ 

This  lemma  reveals  that  the  jump  operator  we  defined  in  (146)  is  in  fact  a  third-order 
perturbation  to  the  classical  energy  jump.  Using  this  fact,  we  can  prove  the  second-order 
accuracy  of  our  numerical  scheme. 


3.7.3  Numerical  dissipation  and  accuracy 

With  the  above  five  lemmas,  we  are  ready  to  state  and  prove  the  main  results  of  the  fully 
discrete  scheme  (152)-(155). 

Theorem  6.  The  solutions  of  the  fully  discrete  scheme  (152) -(155)  satisfy 


H{Phn^ehn+l)-H{pl9hn) 


<0. 


ph .  i  r 

ah 
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A  tr, 


I  dV* 


+  V 


k\W9 


9h  .  ’ n+f 
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T  ,  1  :  VU  .  idV-yr 


n+5  dV. 


eU)' 


lPnf  d3"loc 


9h  ,  Atn  24  <9p3 

n+  4  r 


(Phn+^0hn+l)dvx  + 


3*714 


d3H 


"+3 

(179) 


Proof.  Taking  =  Uh,  in  (152),  wft  =  uh  ,i/9h  .  in  (153),  uA  =  —  1  /  ,  in  (154), 

Tl~\~  2  TL~\~  2  2  ^ 

Wq  =  [p^]/Atn  in  (155),  combing  the  three  equations,  and  following  the  proof  of  Theorem 
4,  one  can  show  that 


bm(Ci;y: 


Uh  ! 

T2  U(  n+2  .Vh  ^ 
13  '  i  In+ 1) 


9h  1 

n+ } 


+  D 
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n+l) 


A  (  Hr  n. 

lAf 
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[p; 
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According  to  Lemma  6,  the  above  relation  can  be  reorganized  as 
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The  last  inequality  is  due  to  Lemmas  3  and  4.  □ 

Remark  19.  This  theorem  implies  that  the  fully  discrete  solutions  respect  the  second  law 
of  thermodynamics.  The  amount  of  dissipation  in  (179)  consists  of  two  parts:  the  physi¬ 
cal  dissipation  and  the  numerical  dissipation.  From  our  analysis,  the  numerical  dissipation 
exclusively  comes  from  the  temporal  scheme,  and  it  will  vanish  if  the  time  step  approaches 
zero. 

Theorem  7.  The  local  truncation  error  in  time  0(f)  =  (@p(t);  ©„(t);  0^(t))T  can  be 
bounded  by  |@(tn)|  <  A'Af2l5  for  all  tn  G  [0, T\,  where  K  is  a  constant  independent  of 
A tn  and  15  =  (1;  1;  1;  1;  1)T . 

Proof.  We  start  by  considering  the  mid-point  rule  applied  to  the  semi-discrete  formulation 
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(138)- (141).  The  fully  discrete  scheme  reads 


Rm  (wh-Yh  )  ■=  (wh 

?  *  n+1/  *  l  ^1  ?  a  » 


V<p!LiU*  x)  =0, 


B^(wh;Yfm)  :=  w\ 


IdXl 
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(v«4.  («A,  +  YelVpJ+.|2  +  /K+.'2 


-Wed+d»UV 
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The  local  truncation  errors  associated  with  the  mid-point  rule  can  be  obtained  by  replacing 
the  time  discrete  solutions  with  the  corresponding  exact  time  continuous  solution: 

B:>;;Au.).u‘(w.),»‘(w.).vi,)  =  K,e—)n, 

Btd(w'‘;/(i„+I),u*((„+1),#'1((„+1),t^J  =  (w‘,@J“)n, 

B^sAw.)V(w).«*(w).f2«)  =  W.es“)„. 


<,  KL 


0h(t„+i) 


|uh(t  ,  l)p 

^c(/(Wi).^(Wi)) - 3^- 
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We  9h{tr 


;Vph(tri 
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Assuming  sufficient  smoothness  for  the  time  continuous  solutions,  one  can  show  that 


e7d  =  o(At2n),  ©r  =  o(Ati) i3,  ©r  =  o(a<), 


wherein  I3  =  (1;  1;  1)T.  Now  replacing  the  time  discrete  solutions  with  corresponding  time 
continuous  solutions  in  the  fully  discrete  formulation  (152)-(154): 

BM(whl]p\tn+l)y{tn+l)1eh{tn+1)1v^)  =  «©p)n, 

Bu (wh;  ph(tn+i),uh(tn+1),dh(tn+1),v£+i)  =  (w\0u)n, 
BE(wh5;ph(tn+1),u\tn+1),e\tn+1),V^1_)  =  («,*  eB)n, 


w%,  Vh  1 

6’  n+2/n 


w, 
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0h(tr 


2  (DocO (in),  O'  (tn+i))  +  hoc(PdWl))  0 h(tn+i)) 


(Pn+1  ~  Pnf  d2"lc 


n+1  Hn!  loc(ph(tn),eh(tn+h)) 
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+  I  w6> 


1  Uft(fn)  •  Uh(tn+i) 


n+|. 


We^(fn+i 

Taylor  expansions  lead  to 
1 


=  0 


2  (I/ioc(pA(in),0A(i„+i))  +  =  "ioc(ph(tn+i),0h(tn+i ))  +  O(A^), 
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Qp2  vr  V-'O)"  V^n+i; 

Uh(tn)  •  uh(in+1)  =  |uh(tn+i)|2  +  0(At2). 


The  above  results  lead  to 


<KL)  -«vi,)  =o(A«2). 


n+2/0 


Due  to  Lemma  7,  one  has 

Pn+l£(Pn+l,Un+l>^+l) 


Pn^(Pn5  <  0%)  [phnE(phn,Uhn,6hn)] 
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=  0{Atl 


Combing  the  above  results  gives  us 


(wf,ej„  =  K,e;“)n  =  Ci(A(j), 

(w‘,eu)n  =  (w'‘,@J“)n  +  0(Ai*)l3  =  0(A£)13, 

«.  eE)n  =  «,  es“)n  +  o(Ai„)  =  o(A«y. 


This  completes  the  proof.  □ 

Remark  20.  According  to  the  proof  of  Theorem  7,  we  can  see  that  the  fully  discrete  scheme 
(152)-(155)  is  a  second- order  perturbation  of  the  mid-point  scheme.  This  perturbation  guar¬ 
antees  the  entropy  dissipation  (179). 


4  Benchmark  problems 

In  this  section,  we  use  a  suite  of  benchmark  problems  to  verify  the  theoretical  estimates  we 
made  in  Section  3. 


Table  1:  One-dimensional  manufactured  solution  for  the  thermal  Navier- Stokes- Korteweg 
equations:  Temporal  convergence  rates  at  t  —  0.5. 

Temporal  errors  in  L2  norm  with  polynomial  degree  k  =  2 


At 

1.0  x  10”1 

IIU  -  y  IU>(n) 

8.00  x  10~4 

order 

- 

IIU  -  UtHym) 

1.50  x  10~2 

order 

- 

IIU  -  Uiiw 

6.23  x  10~3 

order 

- 

5.0  x  10~2  1.0  x  10“2 

1.99  x  10~4  7.96  x  10”6 

2.01  2.00 

3.77  x  10”3  1.49  x  10~4 

1.99  2.01 

1.52  x  10”3  5.95  x  10"5 
2.04  2.01 


5.0  x  10~3  1.0  x  10~3 

1.99  x  10“6  7.93  x  10~8 

2.00  2.00 

3.72  x  10“5  1.48  x  10“6 

2.00  2.00 

1.49  x  10~5  5.94  x  10~7 

2.00  2.00 


Temporal  errors  in  H 1  semi-norm  with  polynomial  degree  k  =  3 


At 

1.0  x  10”1  5.0  x  10~2  1.0  x  10~2  5.0  x  10~3  1.0  x  10~3 

m-nV(Q) 

5.03  x  10~3  1.25  x  10~3  5.00  x  10~5  1.25  x  10~5  5.00  x  10“7 

order 

2.01  2.00  2.00  2.00 

\Y2-Y2h\HHn) 

9.59  x  10~2  2.58  x  10~2  9.84  x  10~4  2.46  x  10~4  9.84  x  10"6 

order 

1.89  2.03  2.00  2.00 

IU  -  UVim 

3.04  x  10-2  7.70  x  10~3  2.78  x  10~4  6.95  x  10“5  2.78  x  10"6 

order 

1.98  2.06  2.00  2.00 
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4.1  Manufactured  solutions 

As  first  examples,  we  construct  one- dimensional  manufactured  solutions  for  the  Navier- 
Stokes-Korteweg  equations  to  corroborate  the  time  accuracy  estimate  given  in  Theorem  7. 
The  computations  are  restricted  to  U  =  (0, 1);  the  exact  density,  velocity,  and  temperature 
for  this  problem  are  chosen  as 

p(x,t )  =  0.5  +  0.1  sin(7rt)  cos(27nr), 
u(x,t)  =  sin(7rt)  cos(27nc), 

6(x,t )  =  0.85  +  0.1  sin(7rt)  sin(47nr). 

The  forcing  terms  for  the  balance  equations  are  obtained  by  substituting  the  above  exact 
solutions  into  the  original  strong-form  problem  (88)-(90).  The  Y  variables  can  be  obtained 
as 


Y\  =  p  =  0.5  +  0.1  sin(7rf)  cos(27nr), 
u  sin(7rt)  cos(27ra;) 

2  6  0.85  +  0.1  sin(7rt)  sin(47rx)  ’ 

Y3  =  --  = - - - . 

6  0.85  +  0.1  sin(7rt)  sin(47rx) 

The  dimensionless  numbers  for  this  verification  problem  are  fixed  to  be  Re  =  1.0,  We  = 
1.0,  and  7  =  1.333;  the  dimensionless  thermal  conductivity  is  taken  as  n  =  1.0.  Periodic 
boundary  conditions  are  enforced  for  all  variables.  The  problem  is  computed  with  spatial 
mesh  size  Ax  =  1.0  x  10~3  for  polynomial  degrees  k  =  2  and  3.  The  time  step  sizes  are  taken 
as  1.0  x  10-1,  5.0  x  10-2,  1.0  x  10“2,  5.0  x  10-3,  and  5.0  x  10-4.  In  Table  1,  the  errors  in 
L2-norm  for  the  quadratic  NURBS  solutions  and  the  errors  in  H1  semi-norm  for  the  cubic 
NURBS  solutions  are  summarized,  ft  can  be  observed  that  the  temporal  errors  converge  like 
0(At2)  in  both  cases.  This  confirms  the  time  accuracy  estimate  given  in  Theorem  7. 

4.2  Coalescence  of  two  vapor  bubbles 

In  this  example,  we  consider  a  one- dimensional  problem  with  zero  sources  (i.e.,  b  =  0  and 
r  =  0)  to  verify  the  entropy  dissipation  estimate  given  in  Theorem  6.  The  computational 
domain  is  set  to  be  f!  =  (0,1).  The  initial  conditions  consist  of  two  static  vapor  bubbles 
with  centers  at  points  C\  =  0.39  and  C2  =  0.61.  The  radii  of  the  bubbles  are  set  to  be 
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R\  —  R'2  —  0.1.  The  initial  density  is  given  by  the  following  hyperbolic  tangent  function. 


p0(x)  =  0.1  +  0.25 


tanh 


‘il(x)~-Rlv/Wef  +  tanh  (  Mx)  ~  R 


2vW) 


dj(x)  =  |x  —  Ci|,  for  i  —  1,2. 


The  initial  velocity  is  zero  and  the  initial  temperature  is  90  =  0.95.  Periodic  boundary 
conditions  are  applied  for  all  variables.  The  dimensionless  numbers  are  taken  as  Re  = 
4.0  x  102,  We  =  1.6  x  10,  and  7  =  1.333;  the  dimensionless  thermal  conductivity  is  taken 
as  k  =  1.0.  The  spatial  mesh  consists  of  104  quadratic  NURBS  functions.  The  problem  is 
integrated  up  to  T  =  10.0  with  time  step  sizes  At  =  1.0  x  10~2,  5.0  x  10-3,  2.0  x  10~3,  and 
1.0  x  10~5. 


Figure  6:  Illustration  of  the  thermal  bubble  dynamics.  pVfl  and  pm  are  the  initial  vapor  and 
liquid  densities;  pa  and  ps  are  the  Maxwell  state  at  the  initial  temperature;  p„,T  and  pir 
are  the  vapor  and  liquid  densities  at  time  T  =  10.0. 


The  two  vapor  bubbles  will  merge  together  to  minimize  the  surface  energy.  At  the 
temperature  9  =  0.95,  the  energetically  stable  liquid  and  vapor  densities  are  0.487  and  0.193 
respectively;  the  initial  vapor  and  liquid  densities  are  pVt0  =  0.1  and  0  =  0.6.  Hence,  the 
vapor  phase  will  become  denser  and  the  liquid  phase  will  become  lighter  to  minimize  the 
free  energy.  In  the  meantime,  the  phase  transition  is  accompanied  with  latent  heat  release 


61 


Figure  7:  Coalescence  of  two  bubbles  for  the  one-dimensional  thermal  Navier-Stokes- 
Korteweg  equations:  Density  profiles  at  (a)  t  =  0.0,  (b)  t  =  0.1,  (c)  t  =  0.5,  (d)  t  =  1.0,  (e) 
t  =  2.0,  and  (f)  t  =  10.0. 
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(f) 


Figure  8:  Coalescence  of  two  bubbles  for  the  one-dimensional  thermal  Navier-Stokes- 
Korteweg  equations:  Temperature  profiles  at  (a)  t  =  0.0,  (b)  t  =  0.1,  (c)  t  =  0.5,  (d) 
t  =  1.0,  (e)  t  =  2.0,  and  (f)  t  =  10.0. 
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Figure  9:  Coalescence  of  two  bubbles  for  the  one-dimensional  thermal  Navier-Stokes- 
Korteweg  equations:  Evolution  of  the  discrete  entropy,  (a)  Global  view;  (b)  Detailed  view 
in  the  vicinity  of  t  —  2.49. 
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and  absorption,  which  will  change  the  local  temperature  distribution.  The  shape  of  the  free 
energy  and  the  Maxwell  states  vary  accordingly.  This  coupled  process  will  eventually  reach 
an  equilibrium  state.  This  dynamic  process  is  illustrated  in  a  density-temperature  phase 
diagram  in  Figure  6.  In  Figures  7  and  8,  snapshots  of  the  density  and  the  temperature  are 
depicted  at  times  t  =  0.0,  0.1,  0.5,  1.0,  2.0  and  10.0.  It  is  observed  that  the  initial  interface 
between  the  two  vapor  bubbles  gradually  vanishes,  and  the  vapor  and  liquid  densities  are 
adjusted  to  achieve  the  energy-stable  states.  The  temperature  of  the  system  fluctuates  in 
time.  The  temperature  first  drops  to  about  0.876  at  time  t  =  1.0,  then  it  raises  to  0.898 
uniformly  at  time  t  =  10.0.  The  Maxwell  states  at  9  =  0.898  are  pv  =  0.1403  and  pi  =  0.5546. 
Figure  7  (f)  shows  that  the  density  at  t  =  10.0  is  very  close  to  the  Maxwell  states.  Since  we 
applied  periodic  boundary  conditions,  the  dissipation  relation  is 
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The  discrete  mathematical  entropy  is  plotted  against  time  in  Figure  9  (a).  It  can  be  observed 
that  H(p^v  6^)  monotonically  decreases  with  respect  to  time,  which  confirms  the  theoretical 
estimate  given  in  Theorem  6.  In  Figure  9  (b),  a  detailed  view  of  the  discrete  mathematical 
entropy  in  the  vicinity  of  t  =  2.49  is  provided.  It  can  be  observed  that  the  differences 
between  the  numerical  solutions  and  the  overkill  solution  decrease  with  reductions  of  time 
step  sizes.  To  verify  the  time  accuracy  estimate,  overkill  solutions  were  first  computed  with 
At  =  1.0  x  10  5.  Then  the  computations  were  repeated  with  larger  time  steps  At  =  5.0 x  10"2, 
1.0  x  10~2,  5.0  x  10~3,  1.0  x  10~3  and  5.0  x  10-4.  The  errors  at  time  t  =  1.0  are  listed  in 
Table  2.  It  can  be  observed  that  the  numerical  solutions  converge  optimally  in  time  to  the 
overkill  solutions.  This  again  corroborates  the  theoretical  estimates  given  in  Theorem  7. 


5  Applications 

In  this  section,  we  investigate  the  van  der  Waals  fluid  model  by  performing  simulations  with 
the  numerical  algorithm  developed  in  Section  3. 
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Tabic  2:  Coalescence  of  two  bubbles  for  the  one- dimensional  thermal  Navier-Stokes-Korteweg 
equations:  Temporal  errors  in  L2-norm  at  time  t  =  1.0. 


At 

5.0  x  10"2 

iiu  - 

8.40  x  10"5 

order 

- 

lin  -  nii-m 

2.69  x  10-4 

order 

- 

lin  - 

2.70  x  10~5 

order 

- 

1.0  x  10"2  5.0  x  10"3 
4.02  x  10"6  1.02  x  10 
1.94  1.98 

6.57  x  10~6  1.64  x  10 
2.02  2.00 

1.50  x  10~6  3.83  x  10 
1.92  1.97 


1.0  x  10"3  5.0  x  10"4 
4.14  x  10~8  1.04  x  10~8 
1.99  1.99 

6  6.54  x  10~8  1.64  x  10“8 

2.00  2.00 

7  1.56  x  10"8  3.93  x  10~9 

1.99  1.99 
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5.1  Evaporation  and  condensation 

In  this  example,  we  numerically  investigate  the  dynamics  of  a  single  vapor  bubble  in  the  pres¬ 
ence  of  temperature  increase  or  decrease  on  the  boundary.  In  this  study,  the  computational 
domain  is  restricted  to  a  unit  square  =  (0,  l)2.  The  center  of  the  vapor  bubble  is  located 
at  the  center  of  the  domain,  i.e.,  Cj  =  (0.5,  0.5);  the  radius  of  the  bubble  is  R\  =  0.25.  A 
hyperbolic  tangent  function  is  utilized  to  give  the  initial  density  profile: 

p0 (x)  =  0.3545  +  0.2479  tanh  jlyWj  , 

di(x)  =  |  x  —  C\  | . 

The  initial  velocity  is  set  to  be  zero.  The  initial  temperature  is  given  by 

90(x)  =  0.85,  if  x  e  £2, 

90(x)  =  9bc ,  if  x  G  dil. 

The  boundary  conditions  for  this  problem  are 

Vp  •  n  =  0,  on  dtt  x  (0,  T), 

u  =  0,  on  x  (0,  T), 

9  =  9bc,  on  x  (0 ,T). 

It  is  known  that  the  hyperbolic  tangent  function  is  only  an  approximation  of  the  steady 
state  solution.  In  the  function  (181),  the  liquid  density  is  0.6024  and  the  vapor  density  is 
0.1066,  which  are  very  close  to  the  Maxwell  states  at  temperature  9  =  0.85.  Hence,  there  will 
be  a  low-intensity  velocity  held  generated  near  the  interfacial  region  to  adjust  the  interface 
profile.  A  snapshot  of  the  velocity  streamlines  at  time  t  =  15.0  is  depicted  in  Figure  10.  If 
9bc  7^  0.85,  as  time  evolves,  thermal  diffusion  will  drive  the  temperature  inside  to  9bc.  This 
change  of  temperature  directly  leads  to  the  change  of  the  Maxwell  states,  which  is  observed 
as  condensation  or  evaporation  of  the  bubble.  If  9bc  >  0.85,  the  liquid-vapor  density  ratio 
will  become  smaller;  if  9bc  <  0.85,  the  liquid-vapor  density  ratio  will  become  larger. 

In  Table  3,  the  Maxwell  states  at  different  temperatures  are  listed.  With  these  values, 
the  radius  of  the  vapor  bubble  at  the  new  stable  configuration  can  be  estimated  by  using 
the  mass  conservation  relation.  Assuming  the  interfacial  region  has  measure  zero,  then  the 


(181) 

(182) 
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total  mass  in  17  is 


0.1066  x  0.2527r  +  0.6024  x  (1.0  -  0.252t r)  =  0.6024  -  0.031tt.  (183) 

If  the  Maxwell-state  liquid  and  vapor  densities  at  the  temperature  9  are  denoted  as  p9  and 
p®,  the  new  radius  of  the  vapor  bubble  Rst  can  be  determined  by  the  mass  conservation 
relation 


p9  x  R2st 7T  +  pf  x  (1.0  -  R2stn)  =  0.6024  -  0.031tt,  (184) 

if  p9  and  pev  satisfy  p9v  <  0.6024  —  0.0317T  <  pf .  If  pf  <  0.6024  —  0.0317T,  the  steady  state  will 
be  a  uniform  liquid  state  with  density  0.6024  —  0.0317T;  if  0.6024  —  0.0317T  <  p9n  the  steady 
state  will  be  a  uniform  vapor  state  with  density  0.6024  —  0.0317T.  The  solutions  of  Rst  for 
9  =  0.95,  0.90,  0.85,  0.80,  and  0.75  are  listed  in  Table  3. 


6 

0.95 

0.90 

0.85 

0.80 

0.75 

Pv 

0.1930 

0.1419 

0.1066 

0.0799 

0.0591 

Pi 

0.4872 

0.5524 

0.6024 

0.6442 

0.6808 

Rst 

- 

0.1916 

0.2500 

0.2802 

0.3000 

Table  3:  The  liquid  and  vapor  densities  at  the  Maxwell  states  of  the  van  der  Waals  fluid  model 
at  different  temperatures.  The  values  are  rounded  to  four  decimal  places.  Rst  represents 
the  steady-state  vapor  bubble  radius  with  the  given  initial  density  profile  (181)-(182)  in  the 
sharp  interface  limit.  When  9  =  0.95,  a  uniform  liquid  state  with  density  p  =  0.5050  will 
form  at  the  steady  state. 

In  the  numerical  simulations,  the  dimensionless  numbers  are  taken  as  Re  =  1.451  x  103, 
We  =  5.263  x  1  and  7  =  1.333;  the  dimensionless  thermal  conductivity  is  taken  as 
k  =  1.378  x  10-3.  The  external  body  force  b  and  the  heat  source  r  are  fixed  to  be  zero.  The 
spatial  mesh  is  comprised  of  5122  quadratic  NURBS  elements.  The  simulation  is  integrated 
up  to  T  =  200.0  with  time  step  size  At  =  1.0  x  10~3.  In  Figure  11,  the  density  profiles 
at  time  t  =  200.0  are  depicted  for  9i)C  =  0.75,  0.80,  0.85,  0.90,  and  0.95.  The  solution  for 
9bc  =  0.95  at  t  =  200.0  forms  a  uniform  liquid  state.  In  Figure  12,  the  density  fields  are 
sampled  over  a  straight  line  y  =  0.5  and  plotted  as  a  single  variable  function  of  x.  Given 
the  density  data  over  the  straight  line,  we  may  introduce  two  points  0  <  27  <  x2  <  1,  such 
that  p{x  1)  =  p{x 2)  =  ( pv  +  pi)/ 2.  In  other  words,  X\  and  x2  are  located  at  the  centers  of  the 
diffuse  interfaces  over  the  line  y  =  0.5.  Then  the  bubble  radius  of  the  numerical  solution  can 
be  defined  as  Rst  :  =  |7'i  —  ^2 1/2.  In  Table  4,  the  values  of  Rst  and  Rst  are  listed.  It  can  be 


Figure  11:  Density  profiles  of  a  single  bubble  under  different  temperature  boundary  condi¬ 
tions:  (a)  Initial  condition,  (b)  9bc  =  0.75,  (c)  9bc  =  0.80,  (d)  9bc  =  0.85,  (e)  9bc  =  0.90,  and 
(f)  9bc  =  0.95. 
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observed  that  the  numerical  results  and  the  theoretical  estimates  of  the  radii  matched  well. 
From  9bc  =  0.90  to  9bc  =  0.75,  the  radius  of  the  bubble  increases  with  the  decrease  of  the 
boundary  temperature. 


Figure  12:  Plots  of  the  density  fields  sampled  over  the  straight  line  y  =  0.5  for  different 
temperature  boundary  conditions.  The  colored  squares  delimits  the  region  (xi,x2). 


9 

0.90 

0.85 

0.80 

0.75 

Rst 

0.1916 

0.2500 

0.2802 

0.3000 

Rst 

0.1902 

0.2492 

0.2798 

0.2998 

Table  4:  The  theoretical  estimates  of  the  bubble  radius  Rst  and  the  numerical  results  of  the 
bubble  radius  Rst. 

To  quantify  the  width  of  the  diffuse  interface  at  a  specific  temperature,  we  introduce  two 
more  points  0  <  X\  <  x2  <  0.5,  such  that  p(x\)  =  Pi~5%(pi—pv )  and  p(x2)  =  Pv+5%(pi~ Pv)- 
Within  the  region  (ah,  £2),  the  density  field  varies  90%  of  the  total  variation  between  the 
liquid  state  and  the  vapor  state.  With  those  two  points  given,  we  may  define  the  interface 
width  as  W  :=  x2  —  The  widths  W  are  listed  in  Table  5.  The  interface  width  increases 
with  the  increase  of  the  temperature.  So  far,  there  has  been  few  results  on  the  estimate  of  the 
interface  width  for  the  van  der  Waals  fluid.  Invoking  the  least  square  method,  we  can  obtain 
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ebc  0.90  0.85  0.80  0.75 

W  0.0123  0.0100  0.0083  0.0070 


Tabic  5:  The  interface  width  W  of  the  numerical  solutions  at  different  temperatures.  The 
values  are  rounded  to  four  decimal  places. 

polynomial  functions  that  best  fit  the  given  data,  which  may  serve  as  empirical  formulas  for 
the  diffuse-interface  width  estimate.  The  linear  polynomial  that  best  fits  the  data  in  the 
least  square  sense  is  W  =  0.0354 9bc  —  0.0198.  The  corresponding  quadratic  polynomial  is 
W  =  0.1010 6lc  —  0. 13136^  +  0.0487.  Both  polynomials  are  plotted  in  Figure  13. 


Figure  13:  The  values  of  the  interface  width  W  at  different  temperatures  are  plotted  as 
red  hexagrams.  These  values  are  fitted  by  linear  and  quadratic  polynomials.  The  linear 
polynomial  is  plotted  as  a  black  solid  line;  the  quadratic  polynomial  is  plotted  as  a  blue  solid 
line. 


5.2  Interface  motion  induced  by  a  temperature  gradient 

It  has  long  been  observed  that  an  externally  imposed  temperature  gradient  may  induce  the 
motion  of  fluid  interfaces.  For  multicomponent  flows,  this  phenomenon  is  ascribed  to  the 
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imbalance  of  surface  tension,  which  may  be  caused  by  differences  in  temperature  [64],  This 
effect  is  referred  to  as  the  thermocapillary  convection  and  is  critical  in  understanding  many 
complicated  physical  phenomena,  such  as  boiling  [49]  and  welding  [50] .  In  the  seminal  paper 
of  Young  et  al.  [75],  surface  tension  is  modeled  as  a  function  of  the  temperature.  This 
theory  was  coupled  with  multiphase  flow  solvers  to  simulate  the  thermocapillary  motion  [2, 
30,  36,  71].  However,  the  interface  motion  for  a  single-component  fluid  under  a  temperature 
gradient  has  rarely  been  studied.  Recently,  a  mathematical  model  was  constructed  based  on 
the  van  der  Waals  theory  [52,  54],  In  these  works,  the  Korteweg  stress  included  an  extra 
term  involving  the  temperature  gradient,  and  the  interstitial  working  flux  II  was  ignored. 
Hence  the  models  they  considered  are  totally  different  from  the  model  we  derived.  In  this 
section,  we  investigate  the  interface  motion  of  a  single  vapor  bubble  under  a  temperature 
gradient  in  two-  and  three-dimensions  with  zero  gravity.  It  is  worth  emphasizing  that  in  our 
simulations  the  capillarity  coefficient  A  remains  a  constant  independent  of  the  temperature. 

5.2.1  Two-dimensional  bubble  motion  in  a  temperature  gradient 

In  this  example,  the  computational  domain  is  a  two-dimensional  square  12  =  (0,  l)2.  The 
boundary  of  12  is  partitioned  into  three  non-overlapping  subdivisions: 

<912  =  r„  u  rt  u  rfe, 

T„  :=  <912  fl  {{x  G  M2|a:  =  0} U {x  G  M2|a:  =  l}}  , 

Tb  :=  <912  fl  {x  G  M2| y  =  0}  , 

Tt  :=<912n{xGM2|r/  =  1}. 

The  initial  density  is  given  by 

Po(x)  =  0.35  +  0.25  tanh  ^  ^ - VWe 

di(x)  =  |x  —  C\  | , 

wherein  the  center  of  the  static  vapor  bubble  is  C\  =  (0.5,  0.5).  The  initial  velocity  is  fixed 
to  be  zero  and  the  initial  temperature  is 

d0(x)  =  o.85  x  g  12  u  rt,  U  Tb, 
d0(x)  =  0.87  x  G  Tt. 
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The  boundary  conditions  for  this  problem  are 


Vp  •  n  =  0, 
u  =  0, 

9  =  0.85, 

9  =  0.87, 

—  q  •  n  =  0, 


on  <911  x  (0,  T), 
on  <912  x  (0,  T), 
on  Tb  x  (0,  T), 
on  Tt  x  (0,  T), 
on  r„  x  (0,  T). 


The  dimensionless  numbers  and  the  dimensionless  thermal  conductivity  are  chosen  as 

Re  =  1.738  x  104, 

We  =  3.277  x  106, 

7  =  1.333, 
k  =  3.453  x  10'3. 


(a)  (b) 

Figure  14:  The  motion  of  a  single  bubble  in  a  two-dimensional  square:  Initial  conditions  for 
density  (a)  and  temperature  (b). 

The  spatial  discretization  is  comprised  of  10242  quadratic  NURBS.  The  time  integration  is 
performed  with  a  fixed  step  size  At  =  5.0  x  10~4  up  to  the  final  time  T  =  500.0.  In  Figure 
14,  the  initial  density  and  temperature  profiles  are  illustrated.  In  Figures  15-16,  the  density, 
temperature,  and  velocity  fields  are  depicted  at  various  time  steps.  It  is  noted  that  there 
is  a  velocity  field  generated  immediately  after  the  simulation  starts.  The  velocity  drives 
the  vapor  bubble  toward  the  positive  thermal  gradient  direction.  Eventually,  the  vapor 
bubble  attaches  to  the  heated  wall  boundary  and  forms  a  hemispheric  shape,  as  is  shown  in 
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Figure  15:  The  motion  of  a  single  bubble  in  a  two-dimensional  square:  Solutions  at  t  —  50.0 
(left  column)  and  t  =  200.0  (right  column).  The  first  row  depicts  the  density  profiles;  the 
second  row  depicts  the  temperature  profiles;  the  third  row  visualizes  the  velocity  streamlines. 
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Figure  16:  The  motion  of  a  single  bubble  in  a  two-dimensional  square:  Solutions  at  t  —  300.0 
(left  column)  and  t  =  500.0  (right  column).  The  first  row  depicts  the  density  profiles;  the 
second  row  depicts  the  temperature  profiles;  the  third  row  visualizes  the  velocity  streamlines. 
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Figure  16  (b).  It  can  be  observed  that,  in  the  liquid  phase,  there  is  a  temperature  gradient 
generated  between  the  heated  top  boundary  and  the  cooled  bottom  boundary.  Inside  the 
vapor  bubble,  the  temperature  distribution  remains  homogeneous  throughout  the  whole 
process.  The  homogeneous  temperature  inside  the  vapor  bubble  was  analyzed  in  [52]  and 
was  attributed  to  the  latent  heat  diffusion. 

5.2.2  Three-dimensional  bubble  motion  in  a  temperature  gradient 

As  a  second  example,  a  three-dimensional  numerical  simulation  is  performed.  The  compu¬ 
tation  domain  is  hi  =  (0,  0.5)  x  (0,  0.5)  x  (0, 1).  The  boundary  of  is  partitioned  into  three 
non-overlapping  subdivisions: 

on  =  rturbu  rvr 

Tb  :=  <912  n  {x  G  R3\z  =  o} , 

Tt  :=  <9f2  Pi  {x  G  R3\z  =  l}  , 

r„  :=  <9f2  fl  {{x  G  M3|a;  =  0} U {x  G  M3|a;  =  0.5} 

U  {x  G  M3| y  =  0}  U  {x  G  M3| y  =  0.5}}  . 

The  center  of  the  vapor  bubble  is  initially  located  at  C\  =  (0.25,0.25,0.3),  and  the  bubble 
radius  is  0.2.  The  initial  density  and  velocity  are 

Po(x)  =  0.35  +  0.25  tanh  ^  ^ - v^We 

di(x)  =  |  x  —  C\  | , 
u0(x)  =  0. 

The  initial  temperature  is 

&o(x)  =  0.85,  x  G  fl  U  Fy  U  Tb, 

6*0(x)  =  0.87,  x  G  Tt. 
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The  boundary  conditions  for  this  problem  are 


Vp  •  n  =  0, 
u  =  0, 

9  =  0.85, 

9  =  0.87, 

—  q  •  n  =  0, 


on  <9f2  x  (0,  T), 
on  dQ  x  (0,  T), 
on  Tb  x  (0,  T), 
on  Tt  x  (0,  T), 
on  r„  x  (0,  T). 


The  dimensionless  numbers  and  the  dimensionless  thermal  conductivity  are  taken  as 

Re  =  3.570  x  103, 

We  =  1.383  x  105, 

7  =  1.333, 
k  =  1.681  x  10”2. 


The  spatial  mesh  for  this  problem  is  comprised  of  128  x  128  x  256  quadratic  NURBS  elements. 
The  time  integration  is  performed  up  to  the  final  time  T  =  200.0  with  a  fixed  time  step 
size  of  At  =  1.0  x  10-3.  Figures  17-19  present  snapshots  of  density  isosurfaces,  velocity 
streamlines,  and  temperature  contours  at  various  slices.  As  soon  as  the  simulation  starts, 
there  is  a  temperature  gradient  generated  in  the  liquid  phase;  the  temperature  field  inside  the 
vapor  bubble  remains  nearly  homogeneous.  Similarly  to  the  two-dimensional  case,  there  is  a 
velocity  field  generated  instantaneously  after  the  top  boundary  is  heated.  The  initial  vapor 
bubble  is  then  driven  by  the  velocity  toward  the  heated  boundary.  At  about  t  =  160,  the 
vapor  bubble  touches  the  top  heated  boundary.  At  t  =  200,  a  vapor  layer  is  formed,  which 
separates  the  heated  wall  boundary  from  the  bulk  liquid  phase.  The  velocity  magnitude  at 
t  =  200  is  uniformly  small.  The  solutions  shown  in  Figure  19  (c)  and  (d)  can  be  regarded  to 
be  very  close  to  the  steady  state  solutions. 


5.3  Boiling 

Boiling  is  a  thermally  induced  phase  transition  process  in  which  new  liquid-vapor  interfaces 
are  generated  in  a  bulk  liquid  region  [15,  57].  The  new  interfaces  may  form  from  discrete 
cavities  on  heated  surfaces,  which  is  called  nucleate  boiling,  or  from  a  stable  superheated 
vapor  layer,  which  is  referred  to  as  film  boiling.  Nucleate  boiling  is  characterized  by  isolated 
bubble  generation  and  is  the  most  efficient  mode  in  heat  transfer.  If  the  surface  temperature 
increases,  bubbles  on  the  surface  tend  to  move  horizontally  and  merge  together  to  form  a 
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Figure  17:  Three-dimensional  motion  of  a  single  bubble:  (a)  density  isosurface  at  t  —  0.0,  (b) 
temperature  on  various  slices  at  t  —  0.0,  (c)  density  isosurface  and  streamlines  at  t  =  40.0, 
and  (d)  temperature  on  various  slices  at  t  =  40.0. 
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Figure  18:  Three-dimensional  motion  of  a  single  bubble:  (a)  density  isosurface  and  stream¬ 
lines  at  t  —  80.0,  (b)  temperature  on  various  slices  at  t  —  80.0,  (c)  density  isosurface  and 
streamlines  at  t  =  120.0,  and  (d)  temperature  on  various  slices  at  t  =  120.0. 
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Figure  19:  Three-dimensional  motion  of  a  single  bubble:  (a)  density  isosurface  and  stream¬ 
lines  at  t  —  160.0,  (b)  temperature  on  various  slices  at  t  =  160.0,  (c)  density  and  streamlines 
isosurface  at  t  —  200.0,  and  (d)  temperature  on  various  slices  at  t  —  200.0. 
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vapor  layer.  Beyond  a  certain  critical  surface  temperature,  a  stable  vapor  film  may  even¬ 
tually  form  between  the  heated  solid  surface  and  the  bulk  liquid  phase,  and  vapor  bubbles 
detach  from  the  layer  periodically.  Film  boiling  is  quite  dangerous  and  should  be  avoided 
in  most  industrial  facilities  because  of  the  heat  accumulated  in  the  vapor  him.  Boiling  has 
been  extensively  employed  in  energy  conversion  facilities,  such  as  power  generators,  cooling 
systems  for  electronic  devices,  and  petroleum  refineries.  Despite  its  importance  in  industry, 
the  fundamental  mechanism  of  boiling  is  still  not  well  understood,  as  was  acknowledged  by 
physicists  [3,  53,  41]  and  engineers  [16,  38].  To  date,  knowledge  about  boiling  is  mainly 
obtained  by  correlating  experimental  data  to  empirical  formulas.  In  view  of  its  disparity  of 
spatiotemporal  scales  and  the  elusive  nature  of  many  subprocesses,  a  predictive  model  for 
boiling  is  highly  important  for  engineering  designs. 

There  have  been  a  few  but  growing  numerical  studies  of  boiling  in  the  past  years.  Film 
boiling  is  regarded  as  most  amenable  to  modeling,  since  its  governing  mechanism  is  princi¬ 
pally  the  Rayleigh- Taylor  instability.  A  multiphase  solver  that  can  simulate  the  Rayleigh- 
Taylor  instability  should  be  capable  of  simulating  him  boiling.  Existing  numerical  simula¬ 
tions  have  been  carried  out  by  the  level  set  method  [67],  the  front  tracking  method  [38],  and 
the  volume-of-huid  method  [73].  Those  simulations  all  started  with  a  preexisting  perturbed 
hat  interface  as  the  initial  condition.  In  other  words,  none  of  those  methods  captured  the 
him  generation  process.  On  the  other  side,  very  few  simulations  of  nucleate  boiling  have  been 
performed  because  more  physical  mechanisms  are  involved  in  this  phenomenon.  A  credible 
nucleate  boiling  solver  is  expected  to  be  capable  of  describing  the  creation  of  new  interfaces 
near  the  nucleation  sites,  handling  the  Rayleigh- Taylor  instability  and  the  Rayleigh- Benard 
instability,  and  tracking  the  moving  interfaces  of  bubbles  and  free  surfaces.  In  [68] ,  the  au¬ 
thors  have  studied  the  nucleate  boiling  by  specifically  designing  a  model  for  the  region  near 
the  nucleate  cavities.  This  approach  destroys  the  conservation  structure  and  relies  on  em¬ 
pirical  data,  including  the  bubble  release  rate,  the  nucleation  site  density,  etc.  In  this  work, 
we  simulate  boiling  flows  in  two  and  three  dimensions,  using  the  Navier-Stokes-Korteweg 
equations.  To  obtain  successful  boiling  simulations,  there  are  several  additional  modeling 
considerations.  First,  the  transport  parameters  are  chosen  to  be  density  dependent  in  order 
to  differentiate  the  properties  of  the  liquid  and  vapor  phases.  Specifically,  the  dimensionless 
viscosity  coefficient  /2  and  the  dimensionless  thermal  conductivity  n  are  larger  in  the  liquid 
region  than  in  the  vapor  region.  In  our  simulations,  these  two  parameters  are  taken  as 


81 


with  Cb°'d  and  Cbml  being  constants  independent  of  p.  Second,  the  gravity  effect  should  be 
taken  into  account  to  generate  the  buoyant  effect.  The  dimensionless  body  force  b  is  chosen 
as 


b  =  (0;  0;  — 0.025)T  , 

for  the  three-dimensional  case  and 

b  =  (0;  — 0.025)T  , 

for  the  two-dimensional  case.  Third,  the  ninety-degree  contact  angle  boundary  condition  is 
used  for  the  density  variable  and  the  slip  boundary  condition  is  applied  to  the  velocity.  To 
specify  the  boundary  condition  for  the  temperature,  the  boundary  dtt  is  divided  into  three 
non-overlapping  parts: 


dn  =  rt  u  rb  u  r„, 

b  =  {x6  <9fi|n(x)  •  b  <  0}  , 

T(,  =  {x6  <9D|n(x)  •  b  >  0}  , 

T„  =  {x  G  <9fi|n(x)  •  b  =  0}  . 

With  the  above  partition,  the  boundary  condition  for  6  is 

6  =  6hl  on  Tj  x  (0,  T), 

0  =  8C,  on  Yt  x  (0,  T), 

— q  ■  n  =  0,  on  T„  x  (0,  T). 

In  the  numerical  calculations,  the  Dirichlet  boundary  conditions  for  0  on  Tj  and  Tt  should 
be  transformed  to  the  Dirichlet  boundary  conditions  for  the  entropy  variable  I5  as 

*5  =  Yb,h  =  on  T6  x  (0,T), 

Vh, 

Yb  =  y5>c  =  on  b  x  (0,  T). 

Vc 

Throughout,  the  Dirichlet  data  are  chosen  as  6h  =  0.950  and  9C  =  0.775.  In  real  situations, 
the  temperature  on  the  solid  surface  cannot  be  evenly  distributed  due  to  surface  unevenness. 
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This  effect  is  modeled  by  adding  perturbations  to  the  Dirichlet  data: 


Y5,h  =  - 

1 

+  5%(x 

0.950 

Wc  =  - 

1 

0.775 

wherein  Sy5  h  (x)  and  8ys  c  (x)  are  scalar  perturbation  functions  that  mimic  the  uneven  tem¬ 
perature  distribution  on  the  solid  surface.  As  for  the  initial  conditions,  the  initial  density 
and  temperature  are  given  by  hyperbolic  tangent  functions;  the  initial  velocity  is  zero.  The 
detailed  formulations  of  the  initial  conditions  are  given  in  subsequent  sections.  These  con¬ 
ditions  represent  a  static  free  surface,  with  liquid  in  the  bottom  region  and  vapor  in  the  top 
region.  It  is  worth  emphasizing  that  the  initial  liquid  and  vapor  densities  are  uniform  with 
no  perturbations.  In  contrast  to  existing  boiling  models,  there  is  no  artificial  manipulation 
used  to  serve  as  boiling  onset  in  this  model.  The  following  results  will  show  that  the  vapor 
bubble  or  the  vapor  film  may  form  automatically  without  preexisting  nuclei  simply  due  to 
local  temperature  variations.  Another  appealing  property  of  this  model  is  that  nucleate 
boiling  and  film  boiling  can  be  simulated  within  a  unified  framework  by  proper  selection  of 
dimensionless  parameters.  These  features  may  be  credited  to  the  thermodynamically  consis¬ 
tent  nature  of  the  model  and  the  algorithm.  It  is  expected  that  this  new  methodology  may 
lead  to  a  predictive  tool  for  boiling. 

The  rest  of  this  section  is  organized  as  follows.  In  Section  5.3.1,  we  perform  a  mesh  sensi¬ 
tivity  test.  In  Sections  5.3.2  and  5.3.3,  two-dimensional  nucleate  boiling  and  two-dimensional 
film  boiling  are  numerically  studied.  A  three-dimensional  boiling  simulation  is  investigated 
in  Section  5.3.4. 

5.3.1  Two-dimensional  mesh  sensitivity  test 

Before  we  start  simulating  practical  examples,  a  mesh  sensitivity  test  is  performed  to  examine 
the  mesh  independence  of  the  solution.  The  simulation  domain  is  11  =  (0, 1)  x  (0,0.5).  The 
material  parameters  are  chosen  as 


We  =  2.103  x  106, 

7  =  1.333, 

cboii  =  2,298  x  10“4, 

Cboti  =  3,448  x  10-5 
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(a)  800  x  400  quadratic  NURBS  elements  (b)  1024  x  512  quadratic  NURBS  elements 


(c)  2048  x  1024  quadratic  NURBS  elements  (d)  4096  x  2048  quadratic  NURBS  elements 


Figure  20:  Density  profiles  of  the  mesh  sensitivity  test  at  t  —  35.0. 


The  initial  conditions  for  this  problem  are 

n  35 

- — - VWe 

u0(x)  =  0, 

00  (x)  =  0.775. 

The  perturbation  for  the  temperature  on  the  boundary  5ys  h  (x)  and  5y5c  (x)  are 

5ysh(x)  =  5.0  x  10~2  sin(107rx), 
hy5  c(x)  =  5.0  x  10-3  sin(107rx). 

The  problem  is  integrated  up  to  T  =  35.0  with  time  step  size  At  =  1.0  x  10~4.  We  use  four 
different  spatial  meshes:  800x400,  1024x512,  2048x1024,  and  4096x2048  quadratic  NURBS 
elements.  The  density  profiles  at  t  =  35  are  depicted  in  Figure  20.  As  can  be  seen,  the 
density  profiles  are  similar  for  all  four  meshes.  In  the  coarsest  mesh,  the  shape  of  the  second 
bubble  attached  to  the  bottom  (from  left  to  right)  is  significantly  different  from  those  in  the 
finer  meshes.  The  solutions  shown  in  Figure  20  (c)  and  (d)  are  almost  indistinguishable. 
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Therefore,  in  the  following  two-dimensional  calculations,  we  used  2048  x  1024  quadratic 
NURBS  elements  to  save  computation  time. 


5.3.2  Two-dimensional  nucleate  boiling 


(a) 


(b) 

Figure  21:  Initial  conditions  of  the  two-dimensional  boiling  simulation:  (a)  density  and  (b) 
temperature. 

In  this  example,  we  simulate  boiling  flows  in  a  two-dimensional  rectangular  domain  hi  = 
(0, 1)  x  (0,0.5).  The  material  parameters  are  chosen  as 

We  =  8.401  x  106, 

7  =  1.333, 

C™  =  1.150  x  10“4, 

cboU  =  1  725  x  10-5 
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(c) 

Figure  22:  Solutions  of  the  two-dimensional  nucleate  boiling  simulation  at  t 
density,  (b)  temperature,  and  (c)  velocity  streamlines. 


=  1.25:  (a) 


(b) 


Figure  23:  Solutions  of  the  two-dimensional  nucleate  boiling  simulation  at  t  —  18.75:  (a) 
density,  (b)  temperature,  and  (c)  velocity  streamlines. 
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(a) 


Figure  24:  Solutions  of  the  two-dimensional  nucleate  boiling  simulation  at  t  —  31.25:  (a) 
density,  (b)  temperature,  and  (c)  velocity  streamlines. 


Figure  25:  Solutions  of  the  two-dimensional  nucleate  boiling  simulation  at  t  —  62.5:  (a) 
density,  (b)  temperature,  and  (c)  velocity  streamlines. 
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Figure  26:  Solutions  of  the  two-dimensional  nucleate  boiling  simulation  at  t  —  100.0:  (a) 
density,  (b)  temperature,  and  (c)  velocity  streamlines. 
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The  initial  conditions  for  this  problem  are 


Po(x) 

=  0.3660  0.2971  tanh  (^2  °'35v/WeJ  , 

(185) 

Uo(x) 

=  0, 

(186) 

0o  (x) 

=  0.775. 

(187) 

In  Figure  21,  the  initial  conditions  for  density  and  temperature  have  been  illustrated.  The 
perturbation  for  the  temperature  on  the  boundary  8y5  h  (x)  and  <5y5  c  (x)  are  uniform  random 
distributions  and  satisfy 


Sy5Jx)  e  [-5.0  x  10~2,  5.0  x  10”2], 

(%,c(x)  e  [-5-0  x  10~3,  5.0  x  10-3]. 

The  spatial  mesh  consists  of  2048  x  1024  quadratic  NURBS  elements.  The  problem  is 
integrated  up  to  the  final  time  T  =  100.0  with  time  step  fixed  as  At  =  5.0  x  10~4. 

In  Figures  22-26,  snapshots  of  the  density,  temperature,  and  velocity  streamlines  are 
depicted.  In  Figure  22,  it  can  be  observed  that  tiny  vapor  bubbles  are  generated  at  discrete 
sites  of  the  heated  wall  surface  during  the  initial  times.  These  small  bubbles  grow  in  size, 
and  some  bubbles  merge  together  to  form  larger  bubbles,  as  is  shown  in  Figure  23.  The 
increase  of  bubble  size  leads  to  the  increase  of  buoyancy.  Beyond  a  certain  critical  point, 
the  bubbles  get  detached  from  the  bottom  boundary  and  rise  upward.  At  about  t  =  18.75, 
the  first  three  bubbles  get  detached  from  the  bottom.  More  bubbles  are  generated  on  the 
bottom  surface  in  the  mean  time.  Figures  24  and  25  show  the  moments  when  two  bubbles 
are  about  to  reach  the  free  surface.  Interestingly,  from  Figures  25  and  26,  small  droplets  can 
be  observed  as  a  result  of  the  breakage  of  the  liquid  film  when  the  vapor  bubbles  reach  the 
free  surface.  There  are  30  bubbles  formed  in  the  time  interval  of  0  <  t  <  100. 

5.3.3  Two-dimensional  film  boiling 

In  the  third  example,  the  same  two-dimensional  problem  considered  in  the  preceding  section 
is  simulated  again  with  a  different  parameter  Cb°tl .  Here,  the  parameter  Cbml  is  chosen  to  be 
4.600  x  10  4,  which  is  four  times  larger  than  that  of  the  previous  example.  Since  the  fluid 
motion  in  this  example  is  slower,  the  simulation  is  integrated  in  time  up  to  T  =  500.0.  All 
the  other  conditions  are  identical  to  those  of  the  previous  case.  In  Figures  27-31,  snapshots 
of  the  density,  temperature,  and  velocity  streamlines  at  different  time  steps  are  presented. 
Once  the  simulation  starts,  a  thin  vapor  film  is  gradually  generated  at  the  bottom  during  the 
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(a) 


(b) 


Figure  27:  Solutions  of  the  two-dimensional  film  boiling  simulation  at  t  —  100.0:  (a)  density, 
(b)  temperature,  and  (c)  velocity  streamlines. 
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(a) 


(b) 


(c) 


Velocity 

0.02«- 


-0.01 


0- 


Figure  28:  Solutions  of  the  two-dimensional  film  boiling  simulation  at  t  —  175.0:  (a)  density, 
(b)  temperature,  and  (c)  velocity  streamlines. 
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(a) 


(b) 


Figure  29:  Solutions  of  the  two-dimensional  film  boiling  simulation  at  t  —  200.0:  (a)  density, 
(b)  temperature,  and  (c)  velocity  streamlines. 
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(c) 

Figure  30:  Solutions  of  the  two-dimensional  film  boiling  simulation  at  t  —  225.0:  (a)  density, 
(b)  temperature,  and  (c)  velocity  streamlines. 
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Figure  31:  Solutions  of  the  two-dimensional  film  boiling  simulation  at  t  —  500.0:  (a)  density, 
(b)  temperature,  and  (c)  velocity  streamlines. 
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early  stage  of  the  simulation  (see  Figure  27).  As  time  evolves,  the  interface  becomes  unstable 
and  there  are  mushroom-shaped  vapor  bubbles  formed,  as  are  shown  in  Figures  28  and  29. 
From  t  =  200.0  to  t  —  225.0,  the  first  two  vapor  bubbles  pinch  off  from  the  vapor  him  and 
rise  upward  in  ellipsoidal  shapes.  As  the  bubbles  get  released  from  the  vapor  him,  two  stems 
are  left  on  the  vapor  him,  which  serve  as  onsets  of  new  bubbles.  This  process  repeats  itself 
periodically.  Till  the  final  time  t  =  500.0,  there  are  totally  seven  bubbles  detached  from  the 
vapor  him.  The  average  bubble  release  rate  for  this  him  boiling  problem  is  much  less  than 
that  of  the  nucleate  boiling  counterpart. 

5.3.4  Three-dimensional  boiling 

As  the  last  example,  we  simulate  the  Navier-Stokes-Korteweg  equations  in  a  three-dimensional 
domain  =  (0, 1)  x  (0,0.5)  x  (0,0.25).  The  material  properties  are  chosen  as 

We  =  6.533  x  105, 

7  =  1.333, 

C\f  =  1.289  x  10~4, 

Cb°il  =  7.732  x  10-5. 

The  initial  conditions  for  this  three-dimensional  problem  are  similar  to  those  of  the  two- 
dimensional  problem,  except  the  free  surface  is  defined  by  X3  =  0.15: 

Po(x)  =  0.33565  —  0.26675  tanh  ^  — ^ — -\/We 
u0(x)  =  0, 

y5)0(x)  =  —1.2334  —  0.0569 tanh  ^  — — — - VWe 

The  perturbation  for  the  temperature  on  the  boundary  5ys  h  (x)  and  5ys  c  (x)  are  uniform 
random  distributions  and  satisfy 


<$yBih(x)  e  [-5.0  x  io~2, 5.0  x  hr2], 

4yV;(x)  e  [-5.0  x  io~3, 5.0  x  io-3]. 

The  spatial  mesh  consists  of  600  x  300  x  150  quadratic  NURBS  elements.  The  problem  is 
integrated  in  time  up  to  T  =  20.0  with  a  fixed  time  step  size  At  =  2.0  x  10“3. 

Remark  21.  The  initial  condition  for  Y5  is  in  fact  a  hyperbolic  tangent  interpolation  of 
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9  =  0.85  for  x3  <  0.15  and  9  =  0.775  for  x3  >  0.15. 


Velocity 


2e-05 


(a) 


Temperature 


0.83  0.95 


(b) 


Figure  32:  Initial  conditions  of  the  three-dimensional  boiling:  (a)  density  isosurfaces  and  (b) 
temperature  isosurfaces. 


In  Figures  32-37,  snapshots  of  density  isosurfaces,  velocity  streamlines,  and  temperature 
isosurfaces  are  presented  at  times  t  =  0.0,  0.6,  5.0,  11.0,  14.0  and  20.0.  At  the  initial 
stage,  there  is  an  unstable  vapor  film  formed  over  the  heated  wall  surface  (see  Figure  33). 
This  film  soon  separates  into  isolated  vapor  bubbles  located  at  random  sites  on  the  bottom 
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(b) 

Figure  33:  Solutions  of  the  three-dimensional  boiling  at  time  t  =  0.6:  (a)  density  isosurfaces 
and  velocity  streamlines;  (b)  temperature  isosurfaces. 
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Figure  34:  Solutions  of  the  three-dimensional  boiling  at  time  t  =  5.0:  (a)  density  isosurfaces 
and  velocity  streamlines;  (b)  temperature  isosurfaces. 
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Figure  35:  Solutions  of  the  three-dimensional  boiling  at  time  t  =  11.0:  (a)  density  isosurfaces 
and  velocity  streamlines;  (b)  temperature  isosurfaces. 
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Figure  36:  Solutions  of  the  three-dimensional  boiling  at  time  t  =  14.0:  (a)  density  isosurfaces 
and  velocity  streamlines;  (b)  temperature  isosurfaces. 
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Figure  37:  Solutions  of  the  three-dimensional  boiling  at  time  t  =  20.0:  (a)  density  isosurfaces 
and  velocity  streamlines;  (b)  temperature  isosurfaces. 
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surface  (see  Figures  34  and  35).  With  the  growth  of  the  bubbles,  the  thermal  energy  is 
conducted  through  the  vapor  region.  Since  the  simulation  domain  is  very  shallow  in  the 
vertical  direction,  these  bubbles  reach  the  free  surface  before  they  get  fully  detached  from 
the  bottom.  When  these  high-temperature  vapor  bubbles  reach  the  cooled  top  surface,  they 
condense  into  liquid  droplets  instantaneously  (see  Figure  36).  At  t  —  20.0,  a  second  round 
of  vapor  bubbles  is  clearly  generated  on  the  bottom  and  the  liquid  droplets  on  the  top 
surface  get  merged  together.  There  is  a  complex  Rayleigh-Benard  mixing  structure  for  the 
temperature  field,  as  is  shown  in  Figure  37. 

6  Conclusions  and  future  work 

In  this  work,  we  presented  a  comprehensive  suite  of  theoretical  and  numerical  methodologies 
for  the  study  of  liquid-vapor  two-phase  flows.  The  contributions  are  elaborated  as  follows. 

A  continuum  mechanic  modeling  framework  for  multiphase  flows  has  been  constructed. 
In  its  derivation,  the  microforce  theory  [27]  is  adopted  with  the  objective  of  accommodating 
non-local  effects.  This  modeling  framework  enjoys  several  appealing  properties.  First,  all 
constitutive  relations  are  represented  in  terms  of  a  thermodynamic  potential.  Therefore,  the 
modeling  work  is  reduced  to  the  design  of  a  proper  form  of  the  thermodynamic  potential. 
Second,  the  framework  automatically  satisfies  the  second  law  of  thermodynamics.  Third, 
some  previously  mysterious  modeling  terms  find  rational  mechanics  explanations  as  a  result 
of  the  Coleman-Noll  procedure.  For  example,  the  “interstitial  working  flux”  [17]  is  the 
power  expenditure  of  the  microstress.  Within  this  framework,  the  Navier- Stokes- Korteweg 
equations  and  the  compressible  Navier-Stokes  equations  are  recovered  by  proper  choices  of 
the  Helmholtz  free  energy  functional. 

A  thermodynamically  consistent  numerical  scheme  for  the  Navier- Stokes- Korteweg  equa¬ 
tions  is  constructed.  For  the  van  der  Waals  fluid  model,  the  definition  of  entropy  variables  is 
generalized  to  the  functional  setting  to  overcome  the  difficulty  induced  by  the  non-convexity 
of  the  entropy  function.  Interestingly,  the  functional  entropy  variables  for  the  van  der  Waals 
fluid  are  formally  identical  to  those  of  the  perfect  gas  model.  The  difference  is  that  the 
entropy  variables  are  not  obtainable  by  merely  an  algebraic  change-of-variables,  but  rather 
they  are  mappings  from  the  conservation  variables  to  their  dual  spaces.  In  the  strong  form,  in 
order  to  obtain  conservation  variables  from  the  entropy  variables,  one  need  to  solve  a  partial 
differential  equation.  An  alternative  statement  of  the  strong  form  problem  is  devised  such 
that  the  equation  for  the  entropy  variable  conjugate  to  density  is  weakly  enforced.  In  doing 
so,  the  weak  formulation  is  guaranteed  to  be  entropy  dissipative,  as  is  shown  in  Theorem  4. 
In  addition  to  the  spatial  discretization,  new  time  integration  schemes  are  developed  based 
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on  a  family  of  new  quadrature  rules  [45].  In  contrast  to  the  traditional  temporal  schemes 
[66],  the  new  time  integration  schemes  do  not  require  convexity  of  the  entropy  function. 
Essentially,  the  new  schemes  can  be  viewed  as  second-order  modifications  to  the  mid-point 
rule.  The  modifications  are  designed  so  that  the  temporal  approximation  is  provably  en¬ 
tropy  dissipative.  The  theoretical  convergence  estimates  have  been  numerically  verified  by 
comparing  solutions  with  manufactured  solutions  and  overkill  solutions. 

The  new  model  and  the  new  algorithm  have  been  applied  to  investigate  a  variety  of  prob¬ 
lems,  including  evaporation,  condensation,  interface  motion  under  a  temperature  gradient, 
and  nucleate  and  him  boiling.  The  advantage  of  the  diffuse-interface  method  is  demon¬ 
strated  by  the  two-  and  three-dimensional  boiling  simulations.  Our  approach  enjoys  several 
desirable  properties.  First,  the  dependency  on  empirical  knowledge  and  assumptions  are 
significantly  reduced.  In  contrast,  existing  boiling  models  rely  heavily  on  empirical  data 
and  sometimes  introduce  artificial  modeling  terms.  Second,  our  approach  provides  a  unified 
modeling  framework  for  both  nucleate  boiling  and  him  boiling.  We  believe  our  methodology 
may  provide  a  predictive  tool  for  a  wide  spectrum  of  the  boiling  phenomena. 

There  are  several  promising  research  directions  for  future  work.  The  Navier-Stokes- 
Korteweg  equations  are  believed  to  be  applicable  to  simulating  cavitating  hows,  the  liquid- 
vapor  phase  transition  induced  by  pressure  variations.  A  potential  challenge  for  such  a 
simulation  is  a  proper  design  of  open  boundary  conditions.  The  van  der  Waals  model  can 
be  further  improved  to  give  more  accurate  material  descriptions.  Recently,  new  models  have 
been  introduced  [65] ,  and  we  anticipate  that  applying  these  new  equation-of-state  may  lead 
to  better  results  in  comparison  with  the  van  der  Waals  model.  On  the  computation  side,  the 
anisotropic  structure  of  the  solutions  makes  mesh-adaptive  isogeometric  analysis  techniques 
[60,  63]  highly  desirable. 
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